J. Phys. Chem. B 2006, 110, 24101-24111
24101
Molecular Dynamics Simulation of the Cybotactic Region in Gas-Expanded Methanol-Carbon Dioxide and Acetone-Carbon Dioxide Mixtures Charu L. Shukla,†,§,| Jason P. Hallett,†,§ Alexander V. Popov,‡,| Rigoberto Hernandez,*,‡,§,| Charles L. Liotta,†,‡,§ and Charles A. Eckert*,†,‡,§ School of Chemical and Biomolecular Engineering, School of Chemistry and Biochemistry, Specialty Separations Center, Center for Computational and Molecular Science and Technology, Georgia Institute of Technology, Atlanta, Georgia 30332 ReceiVed: July 31, 2006
Local solvation and transport effects in gas-expanded liquids (GXLs) are reported based on molecular simulation. GXLs were found to exhibit local density enhancements similar to those seen in supercritical fluids, although less dramatic. This approach was used as an alternative to a multiphase atomistic model for these mixtures by utlilizing experimental results to describe the necessary fixed conditions for a locally (quasi-) stable molecular dynamics model of the (single) GXL phase. The local anisotropic pair correlation function, orientational correlation functions, and diffusion rates are reported for two systems: CO2-expanded methanol and CO2-expanded acetone at 298 K and pressures up to 6 MPa.
I. Introduction Gas-expanded liquids, or GXLs, are a novel and benign class of solvent systems with applications in reactions,1-3 separations,4-6 nanotechnology,7 drug delivery,8,9 and microelectronics.10,11 GXLs are liquid mixtures consisting of an organic solvent combined with a benign gas, such as CO2, in the nearcritical regime. In general, liquid CO2 is a poor solvent, whereas typical organics, such as acetone and methanol, are good solvents. The high compressibility of GXLs enables tunability of physicochemical properties of the liquid solvent, most notably gas solubility, polarity, and dielectric constant. Recent experimental investigations of GXLs have established that replacement of up to 80-90% of an organic solvent with CO2 can enhance reaction efficiency while greatly reducing process waste. Under these conditions, the GXL phase exhibits improved transport properties, dissolves gaseous reactants better, and lead to facile downstream processing. GXLs are a natural offshoot of research involving supercritical fluids. Over the past three decades, supercritical fluids have been used as solvent media for four main reasons: (i) high diffusivities and low viscosities reduce mass transfer limitations, (ii) solvent properties are easily tuned by pressure, (iii) downstream processing and solvent removal are straightforward, and (iv) they provide environmentally benign alternatives to traditional liquid solvents. Supercritical CO2, specifically, has applications in reactions, separations, and microelectronics because of its low surface tension, low reactivity, and high tunability. With no polarity and weak intermolecular forces, however, supercritical carbon dioxide is unable to dissolve many polar and charged species. To some extent, this limitation may be * To whom correspondence should be addressed. E-mail: hernandez@ chemistry.gatech.edu (R.H.);
[email protected] (C.E.). Telephone: (404)894-0594 (R.H.); (404)894-7070 (C.E.). Fax: (404)894-7452 (R.H.); (404)894-9085 (C.E.). † School of Chemical and Biomolecular Engineering. ‡ School of Chemistry and Biochemistry. § Specialty Separations Center. | Center for Computational and Molecular Science and Technology.
mitigated by the addition of a cosolvent under supercritical conditions. Nevertheless, the nearcritical regime, with much better solubilities, would offer a potential alternative to the supercritical regime if only it were sufficiently well characterized that it might be tailored to particular chemical reactions. In this work, we investigate the role of carbon dioxide on altering the cybotactic region, or local environment, in GXLs. The cybotactic region has previously been defined as the region of solute where structure is influenced by solute-solvent interactions.13-16 In Section II, we generalize this definition of the cybotactic region as the region of solvent that is influenced by the presence of a probe. The local composition of a species in such a region is defined as the local number of that species within a specified characteristic radius divided by the total number of molecules within that volume; cf. eq 3 below. The local composition enhancement can thereby be specified as the ratio of the local composition to the bulk composition of a given species. The apparent local composition enhancement of the cosolvent in CO2-expanded methanol is illustrated in the instantaneous snapshot of a MD simulation in Figure 1: while the cosolvent is at low concentration, it is preferentially seen in clusters rather than uniformly distributed within this sample. A more quantitative (and ensemble-averaged) description, provided below, confirms this observation and further suggests that these clusters exhibit a persistence time similar to that in supercritical fluids. As an aid to the reader, Figure 1 also serves the dual service of illustrating the size and relative densities of a typical simulation performed in this work. In supercritical fluids, the cybotactic region seems to exhibit both local composition and local density enhancements due to a combination of solute-solvent interactions (or direct body) and long-range correlations (or indirect n-body interactions) by way of solvent molecules. In this work, we study the cybotactic region in GXLs through a synergy of computation and experiment: experimental data describing the macroscopic (thermodynamic) state are used as inputs to the models, and in turn, the models are simulated to reproduce some of the macroscopic experimental data and elucidate the microscopic structure. This
10.1021/jp0648947 CCC: $33.50 © 2006 American Chemical Society Published on Web 10/31/2006
24102 J. Phys. Chem. B, Vol. 110, No. 47, 2006
Shukla et al. to allow for a solvent molecule to act as the probe, and also to define a cybotactic region independent of the probe, we choose to define Rcut according to the largest distance over which the radial distribution function of the neat solvent displays a significant deviation from the continuum limit. The local density of the solvent within the cybotactic region of a probe is naturally specified by way of this region, i.e.,
Flocal )
N h (Rcut) 4π (R 3 - σexc3) 3 cut
(1)
where N h (Rcut) is the average number of particles of a particular atom type within the solvent sphere of radius Rcut about the probe and given by the following:
Nbulk V
N h (Rcut) ) 4π
Figure 1. Diminishing perspective plot of a snapshot of a simulated trajectory in CO2-expanded methanol at mole fraction CO2 ) 0.88, and P ) 57 bar. Methanol molecules (in blue/shaded) represent the larger structures consisting of three effective atoms: oxygen (Om), protic hydrogen (Hm), and the methyl group (CH3). Because of the diminishing perspective, some of the methanol molecules appear to have only 2 atoms; the third atom is in the linear plane of the molecule. The long pink sticks represent carbon dioxide molecules.
approach is specified in Section III and provides both phenomenological and molecular understanding of the cybotactic region as well as a test for the accuracy of the computational models. In Section IV, the local structure of CO2-expanded methanol and CO2-expanded acetone is illustrated using the radial distribution functions and orientational distribution functions averaged over several molecular dynamics trajectories. The timedependent properties of these systems have also been observed by calculating the local number density autocorrelation functions and diffusion coefficients during these simulations. II. The Cybotactic Region While the use of the term “cybotactic region” has become standard in some disciplines within the literature,13-16 it is not necessarily known to the entire community. Partly for this reason, this section serves to review the use of this term in both the supercritical fluid regime (where it has been used predominantly heretofore) and in the subcritical fluid regime (which is the focus of this article.) In doing so, its definition is generalized so as to accommodate neat fluids and to provide definitions for associated measurables of the cybotactic region. A. Observables in the Cybotactic Region. Differences between the local and bulk compositions have a marked effect on the transport properties and chemical transformations within the cybotactic region. We define the cybotactic region as that spanned by a critical radius, Rcut, about a probe molecule wherein Rcut is specified below. Note that if Rcut were defined according to the distance just beyond the first solvent shell around a spherical probe, then the cybotatic region would be synonymous with the first solvation shell. But this is usually not the case. As summarized in Section B, the cybotactic region has generally been defined in the SCF community as that over which the structure of a neat solvent is displaced by a probe with respect to a particular observable. To extend this definition
∫σR
cut
g(r)r2 dr
(2)
exc
Here Nbulk is the total number of bulk particles, V is the volume of the simulation box and σexc is the excluded radius as defined by the impenetrable volume of the probe molecule. The local number density enhancement is the ratio of the local number density to the bulk number density where Fbulk ) Nbulk/V. The local composition, or xloc, in the cybotactic region may be defined in a similar manner as in ref 17:
xAloc(Rcut) )
N h A(Rcut) N h TOT(Rcut)
(3)
where N h A(Rcut) is precisely N h (Rcut) for a given group A, which can either be a chosen set of molecules or a selected set of atoms whose number density is being measured, and N h TOT(Rcut) is the total number of atoms (from all groups) within Rcut. Although the cybotatic region is the focus of this work, we will also probe the first solvation shell as it is necessarily a subspace of the cybotactic region and is an easily defined indicator of the properties of the larger space. As suggested by Tucker and co-workers,18 the first solvation shell is defined similarly to the cybotactic region, but now Rcut is chosen as the minimum after the first peak in the radial distribution function, g(r), The excluded radius, denoted by σexc, is that radius where g(r) first acquires a nonzero value in the simulation output. B. Cybotactic Region in Supercritical Fluids. The majority of previous studies of the cybotactic region in expanded fluids have focused on the supercritical regime. Both experimental and computational probes of the cybotactic region in supercritical fluids are available. For example, local densities 2-5 times greater than the bulk density were observed through solvatochromic shift measurements,19,20 fluorescence spectroscopy,13-16,21-24 electron paramagnetic resonance spectroscopy,25,26 and UV-vis spectroscopy.27 Evidence of local clustering was observed by Eckert and co-workers,28 who confirmed the existence of large, negative partial molar volumes at solute infinite dilution in attractive mixtures near the critical point. These negative partial molar volumes were also obtained by Chialvo and Cummings29 through integral equation theories (correlation function integrals) of dilute pyrene in supercritical carbon dioxide and dilute neon in supercritical xenon systems. On the basis of partial molar volume data, excess CO2 molecules around naphthalene were calculated by Debenedetti using fluctuation theory.30 Debenedetti also observed local clustering about Xe in supercritical Ne mixtures using molecular dynamics simulations.31 Tucker and co-workers18,32 performed molecular
Cybotactic Region in Gas-Expanded Mixtures dynamics simulations on supercritical, Lennard-Jonesian fluids and quantified both local density enhancements as well as solvent relaxation. C. Cybotactic Region in GXLs. A now-standard approach to altering the local properties of supercritical fluids involves the addition of cosolvents. In the subcritical regime, however, the GXL consists primarily of this putative cosolvent; that is, that which was a cosolvent in supercritical fluids becomes the solvent in GXLs and enables enhanced solubilities. For example, reaction rates for the cis-trans isomerization of 4-4′-disubstituted azobenzenes were shown to decrease with added pressure in both supercritical CO2 and supercritical ethane.12 With the addition of 0.5 mol % cosolvent (methanol, dimethylamine), however, a 15-fold increase in the isomerization rate was observed. The tremendous rate enhancement has been attributed to local composition enhancements of the cosolvent about the solute, which provide an enhanced dielectric environment. For solubility of larger molecules such as organometallics and biomolecules, however, more cosolvent may be needed. Regardless, it should be expected that the cybotactic region in GXLs will exhibit substantially altered solubilities. Previous studies also indicate that the microscopic structure in GXLs may be nonuniform. For example, Subramaniam and co-workers33,34 observed oxygen solubilities greater than those in neat organic solvents. Using this information, the authors ran oxidation reactions of 2,6-di-tert-butyl phenol by Co(salen) catalysts in CO2-expanded dichloromethane and observed 99% conversion. This conversion was found to be significantly higher than in both neat dichloromethane and in scCO2. Furthermore, the turnover frequency in the GXL was forty times higher than in scCO2, while the reaction rates were 10 times faster. III. Computational Methods A. Semiempirical GXL Phase Model. The major challenge in describing the dynamic properties of multiphase systems is the need for a system large enough to accommodate each of the phases and each intervening interface. Such a system is cost prohibitive and has led to the development of multiphase approaches, particularly for vapor-liquid mixtures, that are artificial but utilize a statistically correct algorithm to move particles between phases without the presence of an interface. This approach can provide accurate equilibrium properties and vapor-liquid partitioning behavior, but any such statistical method will artificially affect dynamic behavior, such as kinetic and transport properties. Stubbs and Siepmann17 performed Monte Carlo simulations of aggregation studies in methanol-modified supercritical CO2 at temperature ranges between 303.15 and 324.9 K and 40 MPa and 6.4 mol % methanol. These multiphase calculations allow for the prediction of the partitioning of molecules between the two phases. The good agreement between their prediction and experimental data serves to validate the TrAPPE potentials also employed here for CO2. To obtain estimates of the kinetic properties, we suggest the use of a different artificial constraint. Because we are most interested in the solvent properties of the liquid phase of the vapor-liquid system, we have isolated the liquid-phase region, viz., the GXL stripped of the corresponding vapor phase, and constrained it according to known experimental equilibrium properties. Using model potentials developed in the literature for the components within this effective liquid region, the local solvent and transport properties of a GXL can be obtained without the need to rigorously constrain the system to have the correct vapor-liquid partitioning behavior. Provided significant
J. Phys. Chem. B, Vol. 110, No. 47, 2006 24103 TABLE 1: Solvent Conditions Used for Simulations, T ) 298 Ka solvent
xCO2
P (bar)
FL (kg/m3)
methanol methanol methanol acetone acetone acetone
0.130 0.418 0.884 0.313 0.770 0.920
16 45 57 16 45 57
810.3 846.4 861.8 810.0 860.7 868.5
a Pressures are those reported from the experimental data of Chang.35
and permanent structural changes, of which a phase separation would be an extreme example, do not occur during the course of the simulation, this metastable system should suffice to provide the local static and dynamic behavior in the GXL phase. The quality of the potential models, however, is such that the VLE lines are not perfectly in agreement with experiment, nor have they been performed in all of the conditions of interest to this GXL work. Thus, to sidestep these potential errors, we use the correct solvent density from experiment to obtain the partitioning coefficients in the GXL simulations. Such a choice could lead to a different kind of error in the observations if the potentials were inconsistent with the experimental measurements. The macroscopic properties of this semiempirical GXL model must therefore be validated as done below with respect to macroscopic observables; we chose diffusion coefficients as discussed below. Several studies have shown nonlinear behavior in the properties of GXLs as a function of CO2 mole fraction.19,35 At low CO2 mole fractions (0.10-0.35) the properties of GXLs resemble that of the pure liquid. This region is consequently called the “normal” liquid region in this work. At medium (0.36-0.75) to high (0.8-0.95) CO2 mole fractions, the properties of GXLs resemble those of an expanded liquid in which compressibilities increase. We refer to these regions as the “transition” and “dilated” liquid regions, respectively. In what follows, representative conditions from each of these regimes have been selected and used to conduct detailed and extensive MD simulations. Table 1 presents the respective saturated liquid densities and compositions taken from Chiehming et al.35 selected in the simulations. The simulation system represents a single liquid phase with the density and composition determined by these experiments. B. Force Fields and Ensembles. The mole fractions, densities, and temperature for the acetone-CO2 and methanolCO2 GXL systems investigated here are reported in Table 1. Molecular dynamics (MD) simulations were carried out using the DL_POLY software package.36 Molecules have been modeled as rigid collections of atomic sites with specified fixed charges interacting through pairwise-additive, site-site LennardJones and Coulomb forces:
{( ) ( ) }
uij ) 4ij
σij rij
12
-
σij rij
6
+
qiqj rij
(4)
Carbon dioxide pair interactions have been modeled using the TrAPPE potential.37 The J238 and OPLS-derived potentials39 have been used for methanol and acetone pair interactions, respectively. Both of these potentials are united-atom potentials, treating the methyl groups as one entity. Each of the pairwise potentials specifies a representation for the fixed-point charges, and these are assumed to remain fixed in the heterogeneous pairwise Coulomb interactions. The Lennard-Jones cross-term parameters for MeOH-MeOH and acetone-acetone interac-
24104 J. Phys. Chem. B, Vol. 110, No. 47, 2006
Shukla et al.
tions are handled by the mixing rule specified by the respective potentials. The Lennard-Jones cross-term parameters for MeOHCO2, acetone-CO2, and CO2-CO2 interactions, however, are handled by the Lorentz-Berthelot mixing rule uniformly. As described in the Supporting Information, the use of the Lorentz mixing rule would have altered the parameters by 1% or less. The specific values used for all of the potential parameters are also included in Sections S1 and S2 of the Supporting Information. Regardless of the relative composition, each simulated system box is populated by a total of 1000 molecules. Cubic periodic boundary conditions (PBCs) are used throughout, but the length of the system box is scaled to preserve the specified density as per Table 1. Typical box lengths are on the order of 20-50 Å, and more specifically 40-50 Å for the GXL systems simulated. Because of the PBC, the potential interactions must necessarily be cut off to less than or equal to half the box length. Coulombic interactions are handled by the Ewald summation method with automatic parameter optimization (by DL_POLY). The equations of motion were integrated using the Verlet leapfrog algorithm as implemented by DL_POLY using time steps of 1.0 fs. Structural information (at equilibrium) has been obtained using the NVT ensemble with a Nose-Hoover thermostat in which the relaxation constant is set to 0.3. Initial configurations of the system box have been generated using a random packing of the molecules in the periodic box followed by a 4 ps MD simulation at NVT conditions. Equilibrium averages and time-correlation functions have been obtained as follows: Initial configurations are equilibrated for 50 ps of NVT simulation during which statistics are not collected. This is immediately followed by a 200 ps NVT simulation during which statistics are obtained to determine equilibrium structural information such as g(r) and g(r,θ). An ensemble of 10 configurations sampled from the NVT simulation are propagated for 300 ps under NVE conditions to generate the statistics used to obtain dynamic information, such as timedependent correlation functions and diffusion coefficients. This approach also enabled the determination of error in the diffusion coefficients. IV. Results Local structure in the cybotactic region of CO2-expanded methanol and CO2-expanded acetone has been investigated using the radial distribution function, g(r), and the orientational distribution function g(r,θ). Orientational distribution functions are defined in the following way: after a probe molecule is assigned its local coordinate system (say, r, θ and φ), a number ∆N(r,θ) of atoms of a particular type (say, A) is counted in the finite element with coordinates (r,θ) within the volume
∆V(r,θ) ) 2π
(θ+∆θ)/2 (r+∆r)/2 2 r′ sin θ′ dr′dθ′ ∫(θ-∆θ)/2 ∫(r-∆r)/2
(
) 2π r2∆r +
)( (
3
∆r 12
cos θ -
))
∆θ ∆θ - cos θ + 2 2
)
(
(5) ∆N(r,θ) can be found from the equation
∆N(r,θ) ) [A]
(θ+∆θ)/2 (r+∆r)/2 g(r′,θ′)‚ ∫(θ-∆θ)/2 ∫(r-∆r)/2
2π r′2 sin θ′ dr′ dθ′ (6) which is the definition of g(r,θ) (here [A] is the bulk number
Figure 2. Radial distribution functions for Cg-Cg interactions between carbon dioxide molecules (a) obtained in methanol-carbon dioxide simulations at xCO2 ) 0.130 (solid line), xCO2 ) 0.418 (dashed line), and xCO2 ) 0.884 (dot-dashed line), and (b) in acetone-carbon dioxide simulations at xCO2 ) 0.313 (solid line), xCO2 ) 0.770 (dashed line), and xCO2 ) 0.920 (dot-dashed line). The σ of a Cg atom is 2.8 Å, and therefore the effective σ or σeff is the vdW radius of the Cg-Cg bond, again 2.8 Å. Because the intermolecular potential consists of both Coulombic and LJ interactions, the term “σeff” is used here.
density of atoms A, or equivalently the ratio of the number of those atoms in the entire simulation box to the volume of the box). Given ∆V(r,θ) small enough (∆r is equal to 0.05 Å, ∆θ is 6°), eq 6 becomes
∆N(r,θ) ≈ [A] g(r,θ) ∆V(r,θ)
(7)
This formula is used for calculating the orientational distribution function g(r,θ). The radial distribution function g(r) is defined analogously,
∆N(r) ≈ [A] g(r) ∆V(r)
(8)
but the number ∆N(r) of atoms A is counted within a thin spherical layer
∆V(r) ) 4π
(
)
3
(r+∆r)/2 2 ∆r r′ dr′ ) 4π r2∆r + ∫(r-∆r)/2 12
(9)
Note that g(r) can be obtained from g(r,θ) by integration over θ:
g(r) )
∫0π g(r,θ) sin θ dθ.
A. Structural Information. i. CO2-CO2 Interactions. For simplicity, we characterize the CO2 structure in the vicinity of a chosen CO2 molecule using the relative distance vector between the carbon centers, Cg (carbon in carbon dioxide gas) of the probe and solvent. Figure 3 displays the Cg-Cg radial distribution functions, or g(r), of carbon dioxide in CO2expanded methanol (a) and acetone (b). In the methanol solvent, the Cg-Cg g(r) does not change significantly with addition of CO2 from xCO2 ) 0.130 to xCO2 ) 0.418, remaining at a value near 2.0. At xCO2 ) 0.884, however, the value for g(r) drops to
Cybotactic Region in Gas-Expanded Mixtures
Figure 3. Bulk mole fraction vs local mole fraction of methanol around five different solvatochromic dyes at 40°C: (O) ET-33, (4) 4-nitroanisole, (0) N,N-dimethyl-4-nitroaniline, (]) 4-nitroaniline, (3) 4-nitrophenol. Data extracted from ref 19.
approximately 1.75. In the acetone solvent, the addition of CO2 decreases the first peak in g(r) from 2.1 (xCO2 ) 0.313) to 1.85 (xCO2 ) 0.770) to 1.75 xCO2 ) 0.920), disrupting the CO2-CO2 interactions. In addition, the first g(r) peak in both methanol and acetone solvents are greater than 2.0, indicating significant CO2-CO2 interactions. Figure 4a displays the orientational distribution function, or g(r,θ), for Cg-Cg interactions in CO2-expanded methanol (xCO2 ) 0.130, F ) 810.3 kg/m3). (As illustrated in Figure 4b, θ is defined as the angle between the symmetry axis of the probe CO2 and the vector from the probe Cg atom to a given Cg in the cybotactic region.) Three peaks are clearly observable in g(r,θ). The three calculated peaks are broad, and the peaks are at values of θ equal to 21°, 90°, and 165° at 4.23, 4.03, and 4.23 Å, respectively. The existence of three peaks in Cg-Cg g(r,θ) is consistent with the results obtained by Cipriani et al.40 A slight difference, however, lies in the position of the peaks as they found them at 30°, 90°, and 150° in pure, saturated liquid carbon dioxide at F ) 0.0149 molecules/Å3 and T ) 240 K. The similarity between these two results suggests that the local structure of CO2 in GXLs more closely resembles a condensed liquid than a dissolved gas. These g(r,θ) distributions suggest
J. Phys. Chem. B, Vol. 110, No. 47, 2006 24105 near T-shape like structures between CO2 molecules; a finding also consistent with other simulations in the literature.40-43 The cartoon in Figure 4b shows the most probable local structure between Cg-Cg atoms situated at 21°, 90°, and 165° with respect to a central Cg atom on a carbon dioxide molecule. ii. Methanol-Methanol Interactions. By analogy with the use of Cg above, it is convenient to label the oxygen and the protic hydrogen in a given methanol molecule as Om and Hm, respectively. The intermolecular radial distribution function between the Om probe and nonbonded, protic hydrogen Hm in the cybotactic region around a probe methanol are displayed in Figure 5a. We observe a sharp first peak in Om-Hm g(r) that increases with added CO2 pressure at an r/σ ratio of 1.25 or roughly 2 Å, the length of a hydrogen bond. These g(r) results suggest the manifestation of preferential clustering of methanol molecules. This is consistent with simulations performed by Chatzis and Samios44 and Stubbs and Siepmann.17 The clustering of methanol is illustrated in the arbitrary snapshot of a single trajectory shown in Figure 1. Few, if any, methanol molecules appear outside of the clustered regions. Figure 6a displays the g(r,θ) between oxygen in methanol and the nonbonded, protic hydrogen (Om-Hm interactions) obtained in CO2-expanded methanol simulations: xCO2 ) 0.884, F ) 864.8 kg/m3. (As illustrated in Figure 6b, θ is defined as the angle between the vector formed by connecting Om on the probe methanol to its center of mass and the vector from the Om probe atom to any Hm in the cybotactic region.) Again, a large peak at 150° and 2 Å is observable, suggesting that methanol molecules orient themselves in a very specific manner. The most probable local structure between Om-Hm atoms is shown by the cartoon in Figure 6b. iii. Acetone-Acetone Interactions. The structure in the vicinity of a given aprotic acetone probe molecule is obtained using the relative distance vector between the carbon in acetone and oxygen in acetone (Ca-Oa). As shown in Figure 5b, the radial distribution function of Ca-Oa interactions exhibits two nearly simultaneous peaks at an r/σ of 1.4 and 2.0. This corresponds to radii of 4.69 and 6.7 Å, respectively. These peaks reach a maximum at a g(r) of around 1.4, which is significantly
Figure 4. (a) Orientational distribution function of Cg-Cg interactions in CO2-expanded methanol: xCO2 ) 0.130, F ) 810.3 kg/m3. (b) Most probable orientation of Cg atoms based in Figure 4a.
24106 J. Phys. Chem. B, Vol. 110, No. 47, 2006
Figure 5. (a) Radial distribution function between oxygen in methanol and the nonbonded, protic hydrogen (Om-Hm) from CO2-expanded methanol simulations at xCO2 ) 0.130 (solid line), xCO2 ) 0.418 (dashed), and xCO2 ) 0.884 (dot-dashed) bars. Because the Om-Hm interactions are being specified by LJ and Coulombic interactions, there is no specified σ. In the model, the LJ σ is 3.071 Å and 0 Å for methanol and hydrogen, respectively. Therefore, the effective sigma is the average of the Om and Hm sigmas, or 3.071/2 ) 1. 535 Å. (b) Radial distribution function between the aprotic carbon in acetone and the aprotic oxygen in acetone (Ca-Oa) from CO2-expanded acetone simulations at xCO2 ) 0.313 (solid), xCO2 ) 0.770 (dashed), and xCO2 ) 0.920 (dot-dashed) bars. The σ for Ca and Oa are 3.75 and 2.96 Å, respectively. Therefore, the effective sigma or σeff is 3.35 Å.
smaller than the peaks in the g(r) of Om-Hm interactions (Figure 5a). Furthermore, the first peak in the Ca-Oa interactions at an r/σ ratio of 1.4 (or roughly 4.69 Å) does not change significantly with CO2 addition. The orientational distribution function of Ca-Oa interactions for xCO2 ) 0.313, F ) 810.0 kg/m3 is displayed in Section S3 of the Supporting Information. (Here θ is defined as the angle between the vector formed by connecting the Ca on the probe
Shukla et al. acetone to its center of mass and the vector from the probe Ca to any Ca on an acetone in the cybotactic region.) The orientational distribution function shown in Figure S3a also displays the presence of two simultaneous peaks at 0° and r ) 4.7 Å and r ) 6.7 Å. An additional peak is observed at approximately 100° and 5 Å. Because the two simultaneous peaks are only 2 Å apart and observed at the same angle (0°), it is not physically possible to draw an orientation of these atoms on a two-dimensional surface based on the g(r,θ). Nonetheless, the Ca-Oa orientational distribution function suggests that acetone molecules orient themselves in a specific manner in GXLs. iV. Carbon Dioxide Interactions with SolVent. Radial distribution functions for interactions between carbon dioxide and methanol and carbon dioxide and acetone are presented in Section S4 of the Supporting Information, Figure S4. For CO2expanded methanol simulations, CO2-methanol interactions are represented by the carbon in CO2 and oxygen in methanol (CgOm) g(r). Likewise, we obtain the CO2-acetone structure by analyzing the carbon in CO2 and oxygen in acetone (Cg-Oa) g(r). In CO2-expanded methanol (Figure S4a), the first peak in g(r) at an r/σ of 1.1 Å or approximately 3.2 Å decreases from ∼1 to approximately 0.75 with added CO2. In Figure S4b, a slight increase in the first peak at 3.46 Å from ∼1.42 to ∼1.5, is observed with CO2 addition. This suggests that the correlations between carbon dioxide and acetone are greater than the correlations between carbon dioxide and methanol. Acetone molecules do not cluster nearly as much as methanol molecules, and thus acetone molecules are more available to interact with carbon dioxide. V. Local Number Density. As was rigorously defined in eq 1, the local density (or local number density, in this work) about a probe can be calculated and compared to the bulk density to determine a local (number) density enhancement. Tables 2 and 3 present local number density enhancements for four targeted interactions from CO2-expanded methanol and acetone simulations, respectively. (Here Cg-Cg interactions characterize the local density of Cg atoms around a probe Cg.) While negligible enhancements exist for Cg-Cg interactions, the number density enhancement between Om-Hm interactions increases 6-fold in the dilated liquid region. However, no significant local number
Figure 6. (a) Orientational distribution function of oxygen in methanol and the nonbonded, protic hydrogen (Om-Hm interactions) from CO2expanded methanol: xCO2 ) 0.884, F ) 864.8 kg/m3. (b) Most probable orientation between Om-Hm atoms based in Figure 6a.
Cybotactic Region in Gas-Expanded Mixtures
J. Phys. Chem. B, Vol. 110, No. 47, 2006 24107
TABLE 2: Local Number Density Enhancements (LNDE) for Cg-Cg Interactions and Om-Hm Interactions from CO2-Expanded Methanol Simulations interaction
x (CO2)
Rcut (Å)
Flocal (molecules/Å3)
Fbulk (molecules/Å3)
LNDE (Flocal/Fbulk)
Cg-Cg (normal liquid region) Cg-Cg (dilated liquid region) Om-Hm (normal liquid region) Om-Hm (dilated liquid region)
0.130
5.975
0.0020
0.0019
1.08
0.884
5.975
0.0113
0.0107
1.05
0.130
2.675
0.0141
0.0126
1.11
0.884
2.675
0.0094
0.0014
6.81
TABLE 3: Local Number Density Enhancements (LNDE) for Cg-Cg Interactions and Ca-Oa Interactions from CO2-Expanded Acetone Simulations Fbulk (molecules/Å3)
LNDE (Flocal/Fbulk)
0.0029
0.0028
1.03
5.975
0.0118
0.0106
1.05
0.313
5.425
0.0025
0.0062
0.89
0.920
5.425
0.00947
0.0106
0.89
interaction
x(CO2)
Rcut (Å)
Cg-Cg (normal liquid region) Cg-Cg (dilated liquid region) Ca-Oa (normal liquid region) Ca-Oa (dilated liquid region)
0.313
5.975
0.920
Flocal (molecules/Å3)
density enhancement for Ca-Oa interactions from CO2-expanded acetone simulations is observed. B. Time-Dependent Phenomena. i. Diffusion Coefficients. Self-diffusion coefficients were calculated using the Einstein relation:
D ) lim
1
2 〈[b(t) r - b(0)] r 〉
tf∞ 6t
(10)
Figure 7a plots the simulated self-diffusion coefficients of methanol in CO2-expanded methanol at 323 K alongside those obtained by Aida and Inomata.45 They employ a three-site, flexible model for both methanol and carbon dioxide and incorporate ionization potentials in the mixing rules. Both the J2 potential for methanol used in this work and the flexible potential used in Aida and Inomata’s work were fit to the experimental, saturated liquid density of methanol at ambient temperature and pressure. The subsequent optimized potential
Figure 7. (a) Simulated self-diffusion coefficients of methanol in CO2expanded methanol vs mole fraction carbon dioxide at 323 K: ([) this work, (b) simulated values from Aida and Inomata45 obtained through visual interpolation, and (3) values at 313 K approximated using the Wilke-Chang equation48 with parameters fit from experimental values of Sassiat et al.47 (b) Simulated self-diffusion coefficients of (2) methanol in CO2-expanded methanol and (9) acetone in CO2expanded acetone at 298 K. Self-diffusion coefficients determined experimentally by NMR49 for (0) pure acetone and (4) pure methanol are included. Error bars in the diffusion coefficients were determined by root-mean-square deviations from the average value.
was used to simulate the methanol coexistence curve. Coexistence curves obtained using the methanol potential of Aida and Inomata’s are better correlated to those obtained by experiment.46 As shown in Figure 7a, the simulated self-diffusion coefficients of methanol in CO2-expanded methanol obtained in this work are in reasonable agreement with those obtained by Aida and Inomata at very low and high mole fraction CO2. At less than 50 mol % carbon dioxide, however, the self-diffusion coefficients of Aida and Inomata remain constant with mol % CO2, whereas in the present work, the values steadily increase. For reference, Figure 7a also includes an experimental data point for benzene diffusion in GX-methanol47 at 313 K that has been adjusted to an approximate value for methanol in GX-methanol using the Wilke-Chang equation48 An explanation of the Wilke-Chang ratio is included in Section S5 of the Supporting Information. The simulated self-diffusion coefficients for pure methanol both from this work and Aida and Inomata’s work were less than those determined using Taylor-Aris dispersion. Figure 7b displays the simulated self-diffusion coefficients of (i) methanol in CO2-expanded methanol and (ii) acetone in CO2-expanded acetone. Once again, the results correctly predict an increase in diffusion upon addition of CO2. Peculiarly, the acetone diffusion in the GXL is greater than methanol diffusion at all CO2 concentrations. One expects the converse, as the methanol molecules are smaller both in size and weight compared to the acetone molecules. Again, experimental NMR points obtained by McCall et al.49 are included for reference. While the simulated self-diffusion coefficient for pure methanol at 298 K remains smaller than its experimental point, the simulated pure acetone diffusion is greater. Although our values for diffusion coefficients are not in exact agreement with experimental values, the models reproduce the general trend of enhanced solvent diffusivity with CO2 addition. ii. Local Number Density Autocorrelation Functions. The local number density autocorrelation function is given by
C(t) )
〈δF(0)δF(t)〉 〈δF(0)2〉
δF(t) ) F(t) - Fj
(11)
where F(t) is the local number density at time t, and Fj is the
24108 J. Phys. Chem. B, Vol. 110, No. 47, 2006
Shukla et al. TABLE 4: Instantaneous (τi) and Steady State (τs) Decay Times in the Density Correlation Functions of Figure 8 Are Listed for Several Different Molecular Pairsa interaction Cg-Cg correlations (normal liquid region) Cg-Cg correlations (dilated liquid region) Om-Hm correlations (normal liquid region) Om-Hm correlations ) (dilated liquid region
Figure 8. Local number density autocorrelation functions of: (a) OmHm interactions from CO2-expanded methanol simulations at the dilated liquid region: xCO2 ) 0.884, F ) 864.8 kg/m3; (b) Om-Hm interactions from CO2-expanded methanol simulations at the normal liquid region: xCO2 ) 0.130, F ) 810.3 kg/m3; (c) Cg-Cg interactions from CO2expanded methanol simulations at the dilated liquid region; (d) Cg-Cg interactions from CO2-expanded methanol simulations at the normal liquid region.
ensemble average local number density. The cutoff radius for density determination was taken to be the minimum in the first peak of the radial distribution function, as was employed by Tucker and co-workers.18 Equation 11 can then be discretized to a numerically tractable expression, Ntot
C(t) )
(Ni(Rcut,t) - N h (Rcut))(Ni(Rcut,0) - N h (Rcut)) ∑ i)1 (11)
Ntot
(Ni(Rcut,0) - N h (Rcut))(Ni(Rcut,0) - N h (Rcut)) ∑ i)1 where Ni(Rcut,t) is the local number of neighbor atoms within the cutoff radius at time t, N h (Rcut) is defined in eq 8. Figure 8 displays the local number density autocorrelation functions of both Cg-Cg interactions for a cutoff radius of 5.975 Å in CO2-expanded methanol and Om-Hm interactions for a cutoff radius of 2.675 Å. The local density autocorrelation function (LDAC) is a collective diffusion event, which gives an indication of the local transport in a volume defined by Rcut. Two decay times are obtained from the local density autocorrelation function: the instantaneous (τi) and steady state (τs). These are necessary to distinguish between the initial transient regime and the steady state regime. The former exhibits nonlinear behavior in the mean-squared displacement with respect to time whose decay can be bounded by its initial behavior, whereas the latter regime clearly displays linear behavior. The decay times can be determined from the slope of a semilog plot of the LDAC vs time. Section S6 of the Supporting Information provides the optimal fits to the semilog plots in both the instantaneous regime and steady-state regime of the LDACs shown in Figure 8. Meanwhile, the diffusive escape time for a given molecule in the first solvent shell (or more generally, the cybotactic region) about a probe molecule can be defined as the time it takes for the former molecule to diffuse out of the solvent shell. Assuming diffusive behavior out of a spherical cavity of radius Rcut, this decay time is
Rcut2 tdiff ) 2D
(12)
where D is the diffusion constant. This diffusion constant is
τi (ps)
τs (ps)
diffusive escape time (ps)
11.5
20.4
not available
7.2
14.9
not available
23.5
54.2
12.7
22.2
79.3
5.2
a Diffusive escape time are also included in the two cases wherein diffusion constants were obtained in the simulations (see above)
taken from Figure 7b, and although originally a bulk diffusion constant, it is also assumed to be the same within the first solvent shell. The instantaneous decay, steady-state decay, and diffusive escape times are listed in Table 4. V. Discussion We first address the validity of the models used in this work: the TrAPPE potential37 for CO2 and J238 and OPLSderived potentials39 for methanol and acetone, respectively. It is particularly important in light of our assumption that the models are sufficiently accurate that they would lead to the experimentally observed solvent partitioning if only they were used in a multiphase calculation. Our confidence in the numerical accuracy of this approach is inferred from the following: (1) The simulations never phase separated during our simulations. Thus, the system is at least metastable at the length scales and time scales of the MD simulations explored here. Presumably, if the system were above the spinodal, then it would have phase separated as would be exhibited by complete aggregation of the minor component. (Recall that the system can be in a metastable single-phase “above” the coexistence curve and below the “spinodal” in the sense that the mole fraction of the minority component would be above that at the coexistence curve and below that at the spinodal. In this sense, “above” the spinodal means that the mole fraction of the minority component is higher than that at the spinodal line, and hence the system therein no longer admits a metastable single-phase.) In the vicinity of a mixture’s critical point, the coexistence tie-line is fairly narrow, and hence this requirement is nontrivial though admittedly weak. (2) The calculated bulk diffusion coefficients have the correct order of magnitude and exhibit correct qualitative trends; they are thus in good agreement with the experimental observations of the dense fluid. (3) The results agree qualitatively with the corresponding KamletTaft values (as discussed below), which is an indication that the local averaged structure around the probes is approximately correct. Previous work on supercritical fluids focused on the cybotactic region surrounding an infinitely dilute solute probe. The bulk density of the solvent was therefore equal to the bulk density of the particular supercritical fluid. In GXLs, we have two major components: a solvent and a cosolvent. Therefore, we cannot use the term “bulk density”, as this refers to the bulk density of the total GXL mixture and not individual components. Consequently, we have introduced the terms “bulk number density” and “local number density” in analogy with the composition of the bulk and the cybotactic region, respectively. The presence of carbon dioxide is far more effective in inducing clustering of methanol molecules than acetone mol-
Cybotactic Region in Gas-Expanded Mixtures
J. Phys. Chem. B, Vol. 110, No. 47, 2006 24109
TABLE 5: Kamlet-Tafta Parameters of Methanol, Acetone, and CO2b solvent
R
β
π*
δ (J/cc)1/2
methanol acetone CO2
0.98 0.08 0
0.62 0.51 ∼0
0.59 0.68 ∼0
29.4 19.7
a From refs 52-54. b Parameters for CO2 vary with temperature and pressure, so approximate values are given based on the available literature. See refs 58-60.
ecules. This conclusion can be deduced by analyzing local number density values with CO2 addition. The absolute values of the local number densities of Hm around Om are 0.0141 and 0.0094 molecules/Å3 in the normal and dilated liquid regions, respectively. With a 9-fold decrease in the bulk number density of methanol, we observe relatively little change (∼30%) in the local number density (of Hm around Om). Because these values do not change significantly (despite significant dilution), it is confirmed that methanol forms clusters and remains in clusters in a CO2 environment. However, the absolute values of the local number densities of Oa around Ca are 0.0069 and 0.00104 molecules/Å3 in the normal and dilated liquid regions, respectively. With a 6.7-fold decrease in the bulk number density of acetone, we observe a comparable 6.6-fold decrease in the local number density of (Oa around Ca). In other words, the number of acetone molecules in the first shell is diluted by CO2 molecules in proportion to the mole fraction of CO2. Thus, there is practically no influence due to specific nonideal interactions between the molecules (in contrast to the methanol mixture in which hydrogen bonding strongly distorts the structure). The known intermolecular forces of these molecules50,51 would certainly give similar conclusions, i.e., acetone is aprotic and methanol is well-known to form hydrogen bond oligomers. These results are qualitatively consistent with Kamlet-Taft values. Table 5 presents the Kamlet-Taft52-54 parameters of methanol, acetone, and CO2. The R, or hydrogen bond donating, parameter of acetone (0.08) is low, similar to carbon dioxide. Without hydrogen bonding, there is limited free energy benefit in any sort of self-assembly of the polar species. On the other hand, the R in methanol (0.98) is much greater than CO2. Meanwhile, methanol molecules are not as well solvated in a carbon dioxide environment as they are in hydrogen-bonded clusters. Thus they minimize the system Gibbs free energy by clustering. Through the use of the data obtained by Wyatt et al.19 shown in Figure 2, we have calculated the excess local mole fraction of methanol around several solvatochromic dyes in CO2expanded MeOH at 40°C and various pressures. The local mole fractions were determined by calculating the deviation of the spectral shifts in the GXLs from a linear combination of the values in the pure solvents. These data illustrate the existence of the clustering phenomena we observed in our local structure calculations. The data also illustrate the large variation of the local mole fraction of methanol up to 60% over the bulk mole fraction, with the maximum typically occurring at about 75% CO2 in the mixture. This evidence of local composition enhancement in GXLs is similar to that observed in supercritical CO2-cosolvent systems, but reduced in intensity. Sala et al.55 used IR spectroscopy to study the cybotactic region in CO2-expanded ethanol at 313 K and 10 MPa. They observed an enhancement in solute solubility relative to that in either neat ethanol or pure CO2. Through IR techniques, the authors determined the local structure of ethanol and CO2 molecules about hexamethylenetetramine. Similarly, Kelley, and Lemert56 used a solvatochromic dye, phenol blue, to probe the
cybotactic region in various GXL systems, including CO2expanded methanol and CO2-expanded acetone at 35°C and 55°C and pressures up to 70 bar. These authors found similar effects to those found by Wyatt et al. by utilizing similar assumptions to calculate local compositions.19 The self-diffusion coefficient of pure methanol at 323 K determined in this work (2.85 × 10-9 m2/s) is about 40% of that obtained by Taylor-Aris dispersion techniques57 (6.9 × 10-9 m2/s), which is approximately the uncertainty in Aida and Inomata’s results. This suggests that methanol force fields optimized to fit the saturated liquid density may not be optimal for transport values. However, the significance of our models is that they reproduce accurately the trends we have observed experimentally: that CO2 enhances transport properties (namely diffusion) in the GXL systems investigated. Further, one may expect that the acetone diffusion rate (MW 58) is less than that of methanol (MW 32) because of the lower molecular weight of acetone. However, because of the extent of hydrogen bonding of alcohols in solution,50,51 the “effective” molecular weight of the methanol moieties diffusing may be substantially greater than 58. An alternate explanation, although not explored here, is that the strong hydrogen bonding forces may reduce the effective diffusivity of one methanol molecule. Steady-state time constants obtained from local number density autocorrelation functions show that Om-Hm interactions between methanol molecules persist about 50% longer in the dilated liquid region than in the normal liquid region. Thus CO2 increases the persistence time of methanol interactions, a result which can be explained by the higher R and β Kamlet-Taft parameters of methanol when compared to acetone or CO2. The methanol molecules decrease their Gibbs excess energy in a CO2 environment by interacting with like-like hydrogenbonding forces. As stated in Section IV.B.ii, the local density autocorrelation function gives some indication of local transport events. The longer diffusive decay time seen in the nomal liquid region of methanol suggests that the total number of methanol clusters (or the methanols which exist as clusters) diffuse faster with added CO2. However, the steady-state decay times in the LDAC for the for Om-Hm interactions are smaller in the normal liquid regime than in the diluted liquid regime. This suggests that it takes 25% longer for the methanol to diffuse out of the cybotactic region in the diluted liquid region than in the normal liquid region. The finding is significant, as it shows the persistence time of methanol molecules in the cluster as it relates to the bulk diffusion, a value shown to increase with added CO2 pressure. Three additional conclusions can be drawn from the results of Figure 8 and Table 4: (i) the Om-Hm correlations in both the instantaneous and steady-state regimes persist longer than the Cg-Cg correlations, (ii) the Om-Hm correlations in the dilated liquid region are around 50% longer than those in the normal liquid region for the steady-state regime, and (iii) the Cg-Cg correlations in the dilated liquid region are around 25% less than those in the normal liquid region. Thus, simulations of gas-expanded liquids can be used as a predictive tool. For example, if carbon dioxide is being used as a “miscibility switch” to induce miscibility of two phases,61 MD simulations such as the ones described here can readily be performed to observe the solvation environment of a catalyst as a function of CO2 addition. This provides a feedback to the appropriate CO2 concentration required to optimize the catalysis. VI. Concluding Remarks We simulate CO2-expanded methanol and CO2-expanded acetone at 298 K and several representative densities corre-
24110 J. Phys. Chem. B, Vol. 110, No. 47, 2006 sponding to available experimental data from our lab. Comparison between our simulated diffusion coefficients and experimental values suggests that our GXL models work reasonably well in predicting general trends. Furthermore, the results are qualitatively consistent at the macroscopic level with those found using solvatochromic Kamlet-Taft parameters. The structural results show increased correlations in the radial distribution functions of Om-Hm interactions at around 2 Å with added CO2 concentration. The Om-Hm orientational distribution function also exhibits an aberrant, sharp peak at 150° and 2 Å in the dilated liquid region. These results suggest clustering of protic methanol molecules in CO2-expanded methanol, with the molecules aligning themselves in a specific orientation. Radial distribution functions for Ca-Oa interactions in CO2-expanded acetone, however, do not change with addition of CO2, and thus we do not observe significant clustering of aprotic acetone molecules. In going from the normal to the dilated liquid region, the bulk number densities of acetone (in CO2-expanded acetone) and methanol (in CO2-expanded methanol) decrease 6.7-fold and 9-fold, respectively. A corresponding 6.6-fold decrease in the local number density of Oa around Ca in CO2-expanded acetone is observed, proportional to the dilution process. However, a 9-fold decrease in the bulk number density of methanol does not produce a 9-fold decrease in the local number density. The local number density of methanol molecules remains about the same. Hence, we conclude that carbon dioxide (i) induces methanol clusters in the “normal liquid” region to remain in clusters in the “dilated region” and (ii) directly alters the cybotactic region of acetone molecules in CO2-expanded acetone through a dilution process. Self-diffusion coefficients of acetone in CO2-expanded acetone and methanol in CO2-expanded methanol have both been calculated. Both diffusion coefficients increase with added CO2 concentration; however, the acetone diffusion coefficient is greater than methanol for all CO2 mole fractions simulated. The simulation model accurately reproduces the trends in transport properties, including the diffusion enhancement caused by CO2 addition. Furthermore, time constants obtained through local correlation functions suggest stronger and more persistent interactions between methanol molecules with added CO2 pressure. Identifying structural inhomogeneities in GXL systems may contribute to applications in reacting systems, especially in the design of asymmetric catalysts. For example, the simulation results could help us to understand the solvation shell that surrounds substrates, ligands, and transition states as well as its impact on the reactivity of the critical sites involved in catalysis. Knowing this, we might be able to fine-tune the structure of the ligand or substrates and of their solvation shell to our advantage and favor high enantiomeric excesses. Discerning local structure in GXL microenvironments may also be useful in tuning local solvent properties for size-selective nanoparticle dispersion.62 We feel that this synergistic combination of simulation with laboratory experiments has outstanding potential to advance this understanding, and our future work will focus on determining structural inhomogeneities around a solute molecule. Acknowledgment. This work has been supported by a Department of Energy Basic Energy Sciences grant, DE-FG0204ER15521. C.A.E. acknowledges the support of the J. Erskine Love, Jr., Institute Chair. R.H. is a Goizueta Foundation Junior Professor. The Center for Computational Science and Technol-
Shukla et al. ogy is supported through the National Science Foundation, grant CHE 04-43564. We acknowledge David Bush at Georgia Tech for assistance with equation-of-state models for GXLs. We also acknowledge Jeremy Moix at Georgia Tech for computational assistance in the time-dependent information section. We also acknowledge Mark Maroncelli for providing an early preprint of parallel work in this field63 and for many helpful discussions. Supporting Information Available: Orientation distribution function of Ca-Oa interactions. Radial distribution functions for Cg-Om and Cg-Oa interactions. Explanation of the mixing rules used. Tables of Lennard-Jones parameters and charges used for both CO2-methanol and CO2-acetone simulations. Explanation of the Wilke-Chang ratio used to calculate diffusion coefficients of methanol. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Xie, X.; Brown, J. S.; Joseph, P. J.; Liotta, C. L.; Eckert, C. A. Chem. Commun. 2002, 1156. (2) Jessop, P.; Olmstead, M.; Ablan, C.; Grabenauer, M.; Sheppard, D.; Eckert, C.; Liotta, C. Inorg. Chem. 2002, 41, 3463. (3) Nunes, R. M. D.; Arnaut, L. G.; Solntsev, K. M.; Tolbert, L. M.; Formosito, S. J. J. Am. Chem. Soc. 2005, 127, 11890. (4) Shishikura, A.; Kanamori, K.; Kinbara, H. T. J. Agric. Food Chem. 1994, 42, 1993. (5) Chang, C.; Randolph, A. Biotechnol. Prog. 1991, 7, 275. (6) Sun, Q.; Olesik, S. Anal. Chem. 1999, 71, 2139. (7) Reverchon, E.; Della, G. Powder Technol. 1999, 106, 23. (8) Winters, M. A.; Knutson, B. L.; Debenedetti, P. G.; Sparks, H. G.; Przybycien, T.; Stevenson, C.; Prestelski, S. J. Pharm. Sci. 1996, 85, 586. (9) Bertucco, A.; Lora, M.; Kikic, I. AIChE J. 1998, 44, 2149. (10) Levitin, G.; Myneni, S.; Hess, D. J. Electrochem. Soc. 2004, 151, G380. (11) Myneni, S.; Hess, D. J. Electrochem. Soc. 2003, 150, G744. (12) Dillow, A.; Brown, J.; Liotta, C.; Eckert, C. J. Phys. Chem. A 1998, 102, 7609. (13) Betts, T.; Bright, F. Appl. Spectrosc. 1990, 44, 1203. (14) Zagrobelny, J.; Betts, T.; Bright, F. J. Am. Chem. Soc. 1992, 114, 5249. (15) Zagrobelny, J.; Bright, F. J. Am. Chem. Soc. 1992, 114, 7821. (16) Betts, T.; Zagrobelny, J.; Bright, F. J. Am. Chem. Soc. 1992, 114, 8163. (17) Stubbs, J.; Siepmann, I. J. Chem. Phys. 2004, 121, 1525. (18) Maddox, M.; Goodyear, G.; Tucker, S. J. Phys. Chem. B 2000, 104, 6248. (19) Wyatt, V.; Bush, D.; Lu, J.; Hallett, J.; Liotta, C.; Eckert, C. J. Supercrit. Fluids 2005, 36, 16. (20) Kim, S.; Johnston, K. Ind. Eng. Chem. Res. 1987, 26, 1206. (21) Kajimoto, O.; Futakami, M.; Kobayashi, T.; Yamasaki, K. J. Phys. Chem. 1988, 92, 1347. (22) Knutson, B.; Tomasko, D.; Eckert, C.; Debenedetti, P.; Chialvo, A. ACS Symp. Ser. 1992, 488, 60. (23) Betts, T.; Zagrobelny, J.; Bright, F. ACS Symp. Proc. 1992, 488, 48. (24) Biswas, R.; Lewis, J.; Maroncelli, M. Chem. Phys. Lett. 1999, 310, 485. (25) Carlier, C.; Randolph, T. AIChE J. 1993, 39, 4995. (26) de Grazia, J.; Randolph, T.; O’Brien, J. J. Phys. Chem. A 1998, 102, 1674. (27) Rice, J.; Niemeyer, E.; Dunbar, R.; Bright, F. J. Am. Chem. Soc. 1995, 117, 5832. (28) Eckert, C.; Ziger, D.; Johnston, K.; Kim, S. J. Phys. Chem. 1986, 90, 2738. (29) Chialvo, A.; Cummings, P. AIChE J. 1994, 40, 1558. (30) Debenedetti, P. Chem. Eng. Sci. 1987, 42, 2203 (31) Petsche, I.; Debenedetti, P. J. Chem. Phys. 1989, 91, 7075. (32) Tucker, S.; Maddox, M. J. Phys. Chem. B 1998, 102, 2437. (33) Musie, G.; Wei, M.; Subramaniam, B.; Busch, D. Coord. Chem. ReV. 2001, 219-221, 789. (34) Wei, M.; Musie, G.; Busch, D.; Subramaniam, B. J. Am. Chem. Soc. 2002, 124, 2513. (35) Chiehming, C.; Kou-Lung, C.; Chang-Yih, D. J. Supercrit. Fluids 1998, 12, 223. (36) Smith, W.; Forster, T. R. DL_POLY is a package of molecular simulations routine written by W. Smith and T. R. Forster; The General Council for the Central Laboratory of the Research Councils, Daresbury Laboratory: Daresbury, Nr. Warrington, U.K., 1996.
Cybotactic Region in Gas-Expanded Mixtures (37) Potoff, J.; Siepmann, I. AIChE J. 2001, 47, 1676. (38) Jorgensen, W. J. Phys. Chem. 1986, 90, 1276. (39) Jorgensen, W.; Briggs, J.; Contreras, M. J. Phys. Chem. 1990, 94, 1683. (40) Cipriani, P.; Nardone, M.; Ricci, F. P.; Ricci, M. A. Mol. Phys. 2001, 99, 301. (41) Fedchenia, I.; Schroder, J. J. Chem. Phys. 1997, 106, 7749. (42) Ishli, R.; Okazaki, S.; Okada, I.; Furusaka, M.; Watanabe, N.; Misawa, M.; Fukunaga, T. J. Chem. Phys. 1996, 105, 7011. (43) Zhang, Y.; Yang, J.; Yu, Y. J. Phys. Chem. B 2005, 109, 13375. (44) Chatzis, G.; Samios, J. Chem. Phys. Lett. 2003, 374, 187. (45) Aida, T.; Inomata, H. Mol. Simul. 2004, 30, 407. (46) Honma, T.; Liew, C.; Inomata, H.; Arai, K. J. Phys. Chem. A 2003, 107, 3960. (47) Sassiat, P.; Mourier, P.; Caude, M.; Rossiet, R. Anal. Chem. 1987, 59, 1164. (48) Wilke, C. R.; Chang, P. AIChE J. 1955, 1, 264. (49) McCall, D.; Douglass, D.; Anderson, E. J. Chem. Phys. 1959, 31, 1555. (50) Karachewski, A.; Howell, W.; Eckert, C. AIChE J. 1991, 37, 65. (51) Karachewski, A.; McNeil, M.; Eckert, C. Ind. Eng. Chem. Res. 1989, 28, 315.
J. Phys. Chem. B, Vol. 110, No. 47, 2006 24111 (52) Kamlet, M.; Abbound, J.; Taft, R. J. Am. Chem. Soc. 1977, 99, 6027. (53) Kamlet, M.; Taft, R.; J. Am. Chem. Soc. 1976, 98, 377. (54) Kamlet, M.; Taft, R.; J. Am. Chem. Soc. 1976, 98, 2866. (55) Sala, S.; Ventosa, N.; Tassaing, T.; Cano, M.; Danten, Y.; Besnard, M.; Veciana, J. Chem. Phys. Chem. 2005, 6, 587. (56) Kelley, S.; Lemert, R. AIChE J. 1996, 42, 2047. (57) Frank, M. J. W.; Kuipers, J. A. M.; Swaaij, W. P. M. J. Chem. Eng. Data 1996, 41, 297-302. (58) Sigman, M. E.; Lindley, S.; Leffler, J. E. J. Am. Chem. Soc. 1985, 107, 1471. (59) Yonker, C. R.; Frye, S. L.; Kalkwarf, D. R.; Smith, R. D. J. Phys. Chem. 1986, 90, 3022. (60) Maiwald, M.; Schneider, G. M. Ber. Bunsen-Ges. Phys. Chem. 1986, 90, 960. (61) West, K. N.; Hallett, J. P.; Jones, R. S.; Bush, D.; Liotta, C. L.; Eckert, C. A. Ind. Eng. Chem. Res. 2004, 43, 4827. (62) McLeod, M. C.; Anand, M.; Kitchens, C. L.; Roberts, C. B. Nano Lett. 2005, 5, 461. (63) Li, H.; Maroncelli, M. J. Phys. Chem. B 2006, http://dx.doi.org/ 10.1021/jp064166j.
