J. Phys. Chem. B 1997, 101, 5473-5476
5473
Understanding the Mechanism for the Mass Accommodation of Ethanol by a Water Droplet Ramona S. Taylor, Douglas Ray, and Bruce C. Garrett* EnVironmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, Richland, Washington 99352 ReceiVed: February 21, 1997; In Final Form: May 1, 1997X
The interaction of small gas-phase molecules with the liquid/vapor interface of aqueous droplets and their subsequent accommodation into the bulk of the droplet is an integral part of the chemistry of the troposphere. On the basis of an analysis of a kinetic mechanism for the mass accommodation process, Worsnop, Davidovits, and co-workers (J. Phys. Chem. 1996, 100, 13007) predict that most gas-phase molecules must surmount a substantial free energy of activation before becoming solvated into the aqueous phase. In this Letter, molecular dynamics computer simulations are used in conjunction with statistical mechanical perturbation theory to examine the molecular-level details of this process. Due to the availability of experimental data with which to compare our findings, the interaction of ethanol with a H2O lamella was chosen as a prototypical system for study. The calculated equilibrium free energy surface for transporting ethanol across the liquid/vapor interface and into bulk H2O exhibits a barrier to solvation that is 8.2 kcal/mol smaller than that predicted by the Worsnop/Davidovits model. This discrepancy suggests that nonequilibrium solvation or other kinetic effects may dominate the transport of small molecules across the liquid/vapor interface of water.
The uptake of gas-phase molecules by cloud droplets and aerosols is an essential step in atmospheric processes leading to the removal and transformation of many species. Important processes include the production of acid rain and the removal of anthropogenic species such as the replacements for CFCs. Thus, an understanding of the uptake process, or at least accurate values of the critical parameters, is essential to the development of atmospheric chemistry models with predictive capabilities. The uptake process for nonreactive species consists of diffusion in the gas phase to the liquid/vapor interface, transport across the interface, and diffusion in the aqueous phase away from the interface. A number of elegant experimental techniques have been used to measure the uptake of small gas-phase molecules by aqueous droplets under a variety of conditions.1-3 Transport across the interface begins with a molecule in the gas phase impinging upon the liquid/vapor interface. Because of the large number of molecular modes in the water droplet, many at low frequencies, the water droplet rapidly dissipates the thermal energy of the impinging molecule; thus, it is a good approximation to assume a sticking coefficient of unity. The thermally accommodated molecule can then (i) desorb back into the vapor phase, (ii) remain adsorbed on the surface, or (iii) become solvated in the bulk liquid phase. The existence of stable surface states for typical nondissociating gas-phase species, e.g. ethanol, is expected due to the relative lessening of the surface free energy of the liquid/vapor interface obtained via the addition of these molecules to H2O.4,5 Direct observations of surface states for SO2 and dimethyl sulfoxide on aqueous surfaces6,7 have recently been accomplished by using surface second harmonic generation spectroscopy.8 The probability that a specific gas-phase molecule will be accommodated by a liquid droplet is represented phenomenologically by a mass accommodation coefficient, R
R)
number of molecules entering the liquid phase number of molecular collisions with the surface
(1)
Thus, R equal to 1 indicates that the impinging gas-phase X
Abstract published in AdVance ACS Abstracts, June 15, 1997.
S1089-5647(97)00644-5 CCC: $14.00
molecule will be accommodated by the liquid and not desorb back into the gas. The mass accommodation coefficients for more than 30 small molecules interacting with H2O droplets have been measured.1,3 For all of these molecules, R is less than 0.25 at 273 K.1,3 In addition, the mass accommodation coefficients show a negative temperature dependence, indicating the probable importance of a surface state in the uptake process.1,3 From these and other similarities in the general behavior of the measured mass accommodation coefficients, a model for the mass accommodation process has been presented by Worsnop, Davidovits, and co-workers.2 Within this model, the impinging molecule initially makes contact and, assuming a sticking coefficient of unity, forms a weak interaction with a limited number of H2O molecules at the surface. It can then either desorb back into the gas or be accommodated by the surface and subsequently solvated into the bulk liquid. Consequently, assuming that the concentration of molecules in the bulk liquid is far from saturation, the mass accommodation coefficient, R, is merely a measure of the probability that the molecule will be solvated into the bulk liquid versus desorbed back into the gas phase:1
(
)
-∆Gqsolv exp ksolv -∆(∆Gq) R RT ) ) ) exp 1 - R kdesorb -∆Gdesorb RT exp RT
(
)
(
)
(2)
Here, ksolv and kdesorb are the solvation and desorption rate constants, respectively, ∆Gqsolv and ∆Gdesorb are the activation free energies of solvation and desorption, respectively, and ∆(∆Gq) is the difference in the activation free energies. If the energetics for the process of a molecule desorbing from the liquid back into the gas phase are assumed to be monotonically uphill, then the activation free energy for desorption is equivalent to the free energy for desorption. Given a mass accommodation coefficient of 0.1, a molecule would have an activation free energy difference, ∆(∆Gq), of +1.3 kcal/mol at 273 K. The justification for this barrier is phrased in terms of nucleation theory.2 However, since the present experimental techniques © 1997 American Chemical Society
5474 J. Phys. Chem. B, Vol. 101, No. 28, 1997 typically measure macroscopic quantities, the molecular-level information necessary to validate this theory is not available. In this Letter, we use molecular dynamics (MD) computer simulations9 in conjunction with statistical mechanical perturbation theory10 to explore the equilibrium energetics of the mass accommodation process. As a prototypical system for the mass accommodation process, the ethanol-H2O system was chosen for investigation. As the reaction coordinate, we choose the distance between the centers of mass (COM) of the ethanol and water droplet, such that this coordinate is perpendicular to the interface. Along this reaction coordinate, the effective energies of interaction of these two moieties are averaged over the molecular orientations in the ethanol and water droplet. The local maxima in this free-energy profile can then be used to estimate ∆(∆Gq). This is only an estimate since the calculation assumes that the water molecules stay in equilibrium with the ethanol molecule throughout the accommodation process. MD simulations have previously been utilized to explore the structures and energetics of both the H2O liquid/vapor interface11-13 and the hydrophobic H2O/“organic liquid” interface.14 Statistical mechanical perturbation theory has widely been used to calculate the potential of mean force (PMF) for the association of two entities in a bath of solvent molecules.15 The PMF provides a means to quantitate the effect of the bath molecules on the intermolecular interactions of the two entities along a reaction coordinate. The difference in the PMFs between two successive configurations is equal to the reversible work, W(r), or free energy difference, for this configurational change in the presence of the solvent.15 The relative free energy difference between two configurations can be calculated as follows
W(r) ) A(r0(∆r) - A(r0) ) -kBT ln(〈exp{-β[U(r0(∆r) - U(r0)]}〉r) (3) where the Helmholtz free energy difference is calculated between a reference system, r0, and a system where the reaction coordinate is perturbed from r0 to (r0 ( ∆r). U(r0) and U(r0(∆r) are the potential energies of the reference system and of the perturbed system, respectively. 〈...〉r indicates that the average is calculated relative to the potential energy of the reference state, U(r0). A complete equilibrium free energy curve as a function of the reaction coordinate can be ascertained from a series of such configurations. Equation 3 is defined in a canonical (NVT) ensemble (where the number of molecules, volume, and temperature are held constant) and yields Helmholtz free energies. In the present simulations, the ensemble we are using more closely mimics an isothermal-isobaric (NPT) ensemble (where the number of molecules, pressure, and temperature are held constant), and thus the free energies calculated from our simulations can be directly compared to the experimentally measured Gibb’s free energies. For the calculations reported here, the simulation cell consists of a single ethanol molecule and 525 H2O molecules sandwiched between two (25 Å)3 volumes of vacuum. Periodic boundary conditions are applied in all three directions. All simulations are performed at a temperature of 298 K. The H2O model employed in these calculations is the rigid SPC/E water potential developed by Berendsen and co-workers.16 This potential has been shown to adequately reproduce the surface tension and structural properties of the H2O liquid/vapor interface.11,12 The free energy of solvation for a H2O molecule into bulk SPC/E H2O is calculated to be -7.4 kcal/mol17 relative to the experimental value of -6.3 kcal/mol.18 The intramolecular potential model for ethanol was derived from the AMBER 4.0
Letters force fields.19 The atom-centered charges were calculated via the Gaussian 94 program20 using a Merz-Kollman electrostatic potential fit21 at the MP2/6-31+G level of theory. The intermolecular ethanol-H2O parameters were determined by the Lorentz-Bertholot mixing rules.9 All of the molecular dynamics simulations were run using a modified version of the AMBER 4.0 suite of programs.19 A molecular dynamics step size of 0.001 ps was employed. Heating was accomplished via the Berendsen scheme with a coupling constant of 0.2.22 The bond lengths in both the H2O and ethanol molecule were constrained with the SHAKE algorithm.23 The angles in the ethanol molecule, however, were not constrained. All simulations utilized a cutoff distance of 9 Å. A complete description of the ethanol potential and of the calculation details will be given in full later.17 Initially, the ethanol molecule was placed 2 Å above an equilibrated H2O surface as depicted in Figure 1A. After an equilibration period of 500 ps, the ethanol molecule was found to lie at the H2O surface as shown by the density profile given in Figure 1C. This profile is constructed by dividing the simulation cell into 1 Å thick layers along the z axis and subsequently calculating the ethanol and H2O densities in each layer. Due to the construction of the simulation cell discussed above, the resulting H2O density profile exhibits two interfaces with the plateau region having a density consistent with that of bulk H2O. However, in the free energy calculations discussed below the ethanol molecule traversed across only one of these two interfaces before entering the bulk water; thus, only that half of the density profile which corresponds to the free energy calculations is shown here. The H2O density profile is fit with a hyperbolic tangent function as explained in ref 11. The ethanol density is found at only one surface of the H2O lamella. During the course of these simulations, the ethanol molecule was never found to enter either the gas phase or the bulk aqueous phase. This can be understood from the free energy profile and will be discussed below. Calculations of the orientation of the ethanol molecule with respect to the H2O surface reveal that, as expected, the hydrophobic end of the ethanol molecule points out of the surface. As shown in Figure 2, the angle between the surface normal and the C-CH3 bond is broadly distributed about 63°, which is indicative of the bent geometry of ethanol. The free energy curve for the insertion of an ethanol molecule into bulk H2O has been calculated by using the thermodynamic perturbation approach as described above. The reaction coordinate, z, is as defined above. A series of simulations were performed in which the COM of the ethanol molecule was moved by 0.5 Å along the reaction coordinate toward the bulk water. Each simulation was equilibrated for 50 ps prior to a 100 ps trajectory in which energetic data was collected every 5 fs. The free energy curve constructed from this series of simulations is given in Figure 1B. The free energy differences calculated via the statistical mechanical perturbation method are relative and not absolute. Here, the zero of energy is defined to be when the ethanol molecule is far from the H2O surface. Thus, the free energy curve can be adjusted such that a measure of the absolute free energies of desorption, ∆Gdesorb, and solvation, ∆Gsolv, can be calculated for ethanol in H2O as defined by the arrows in Figure 1B. (∆Gsolv is a measurement of the free energy difference needed to solvate a molecule relative to the gas-phase energy and should not be confused with ∆Gqsolv of eq 2, which is the activation free energy for solvation and is measured relative to the free energy of a molecule residing on the H2O surface.) ∆Gdesorb and ∆Gsolv are found to have values of 7.2 (( 0.1) and -5.8 (( 0.1) kcal/mol, respectively. The error bars are 1 standard deviation of five 20 ps trajectories.
Letters
J. Phys. Chem. B, Vol. 101, No. 28, 1997 5475
Figure 2. Average orientation of the ethanol molecule relative to the surface normal. θ is defined as the angle between the surface normal, nz, and either the C-O bond or the C-CH3 bond. Since an isotropic distribution is flat in cos θ, the relative intensities are plotted as a function of cos θ instead of as a function of θ. The open circles correspond to the orientation of the C-O bond, and the closed circles correspond to the orientation of the C-CH3 bond. The inset depicts the typical orientation for an ethanol molecule on H2O. The color scheme for the ethanol molecule is the same as that of Figure 1.
Figure 1. (A) Depiction of the interaction of an ethanol molecule with a portion of the H2O lamella at 298 K. The large gray and small black spheres represent the O and H molecules in H2O, respectively. The small light gray, medium gray, large dark gray, and small gray spheres represent the alkyl H, the C atoms, the O atom, and the alcoholic H atom of ethanol, respectively. (B) Free energy profile for the insertion of an ethanol molecule into bulk H2O. The definitions of ∆Gdesorb, ∆Gsolv, and ∆(∆Gq) used in this paper are shown by the arrows. The heavy solid line indicates the height of the experimentally measured activation barrier between the surface site and the bulk solvated site calculated within the Worsnop/Davidovits model. (C) Density profile of 1 ethanol and 525 H2O molecules equilibrated at 298 K for 500 ps. The open circles correspond to the H2O density profile, with the solid line corresponding to a fit to this data obtained using the hyperbolic tangent function given in eq 2 of ref 11. The squares correspond to the ethanol density profile. Here the solid line is meant only to guide the eye.
Our calculated value for ∆Gsolv agrees well with the experimental value of -5.05 kcal/mol.18 Consequently, we feel confident that our potential correctly describes the H2O-ethanol energetics. These simulations predict that the subsequent solvation of an ethanol molecule that is adsorbed to an H2O surface into bulk water is endothermic by 1.8 kcal/mol. This barrier should restrict the solvation of the ethanol molecule within the 500 ps lifetime of a typical MD simulation to that of a rare event. From the free energy profile (Figure 1B), ∆(∆Gq) is calculated to be -5.4 (( 0.1) kcal/mol. By using this value in eq 2, the mass accommodation coefficient for ethanol in H2O
at 298 K is calculated to be 1.0. The corresponding experimental value is 0.0092.3 An experimental value for ∆(∆Gq) can be obtained by substituting this value of R into eq 1. The resulting experimental ∆(∆Gq) is +2.8 (( 0.5) kcal/mol. This is 8.2 (( 0.5) kcal/mol larger than our calculated value. The potentials used in these simulations correctly reproduce the energetic and dynamic properties of bulk H2O, the surface tension and structure of the liquid/vapor interface of H2O, the energetics of ethanol solvated in bulk water, the structure of ethanol adsorbed at the liquid/vapor interface of H2O, and the preference for ethanol to adsorb at the H2O surface; thus, this 8.2 kcal/mol difference in ∆(∆Gq) is not a result of an inadequate potential model. Our simulations only measure the equilibrium free energy profile; thus, the differences between the experimental value for ∆(∆Gq) and that of our simulations may be due to nonequilibrium effects. For example, reorganization of the solvent molecules into nonequilibrium configurations that will accommodate the ethanol molecule may influence the kinetics. Simulations aimed at quantifying the nonequilibrium effects on the solvation barrier are underway.24 It is unprecedented, however, that classical theories of nonequilibrium solvation dynamics25 will account for this additional 8.2 kcal/mol barrier, which at room temperature would correspond to a 6 order of magnitude decrease in the reaction rate. It should also be noted that in these simulations, the reaction coordinate is chosen to be the distance in the z direction between the center of mass of the box of H2O and the center of mass of the impinging ethanol molecule. In the Worsnop/Davidovits model, at the minimum of the surface site, the reaction coordinate changes from that which we have chosen to one denoting the degree of nucleation of the solvent around the impinging molecule. Due to the constraints which we impose on our definition of the reaction coordinate, the barrier calculated here is, at worst, an upper bound to the true equilibrium solvation barriersnot a lower bound as would be needed to bring our findings into accord with those of the proposed model for mass accommodation. Consequently, additional experimental and theoretical investigations need to be undertaken to examine why these differences exist. Perhaps then, a true molecular-level understanding of the mass accommodation process will be in hand. Recently, Wilson and Pohorille26 have reported complementary calculations of the uptake of ethanol by a H2O surface at 310 K. These calculations used the TIP4P potential model for
5476 J. Phys. Chem. B, Vol. 101, No. 28, 1997 H2O and an ethanol potential that was derived from the AMBER force fields but that has different values for the atom-centered charges than those used here. These authors calculated a free energy surface for the insertion of ethanol into bulk H2O (using the umbrella sampling method) that is similar in shape to that shown here. It, however, predicts values for ∆Gdesorb, ∆Gsolv, and ∆(∆Gq) of 6.3, -5.2, and -4.9 kcal/mol, respectively, compared to the corresponding values of 7.2, -5.8, and -5.4 kcal/mol found in this work. Unlike the SPC/E model, the TIP4P model is known to underestimate the surface tension of water.27 This may be responsible for the differences in energetics between our calculations and theirs. In addition to these equilibrium solvation calculations, Wilson and Pohorille have performed calculations of the scattering of an ethanol molecule from a water surface in which no appreciable dynamical barrier to solvation was found. They have also performed MD simulations of ethanol entering the bulk H2O that show that the ethanol motion is diffusive in nature. These results are consistent with our major conclusion that the calculated free energy barriers do not describe the experimentally observed mass accommodation coefficients and that other dynamical effects still need to be considered. Acknowledgment. We thank Drs. Liem X. Dang, Greg Schenter, Paul Davidovits, Douglas Worsnop, and Michael Wilson for insightful discussions. This work was performed under the auspices of the Division of Chemical Science, U.S. Department of Energy, under contract DE-AC06-76RLO 1830 with Battelle Memorial Institute, which operates the Pacific Northwest National Laboratory. R.S.T. also acknowledges a postdoctoral fellowship from Associated Western Universities, Inc. under Grant No. DE-AC06-76RLO 1830 with the U.S. Department of Energy. References and Notes (1) Kolb, C. E.; Worsnop, D. R.; Zahniser, M. S.; Davidovits, P.; Hanson, D. R.; Ravishankara, A. R.; Leu, M.-T.; Williams, L. R.; Molina, M. J.; Tolbert, M. R. Current Problems and Progress in Atmospheric Chemistry, In AdVanced Series in Physical Chemistry; Ng, C.-Y., Ed.; World Scientific: Singapore, 1995; Vol. 3, pp 771-875. (2) Nathanson, G. M; Davidovits, P.; Worsnop, D. R.; Kolb, C. E. J. Phys. Chem. 1996, 100, 13007. (3) Jayne, J. I.; Duan, S. X.; Davidovits, P.; Worsnop, D. R.; Zahniser, M. S.; Kolb, C. E. J. Phys. Chem. 1991, 95, 6329. (4) MacRitchie, F. Chemistry at Interfaces; Academic Press, Inc.: New York, 1990, pp 4-23. (5) Dynarowicz, P. Colloid Polym. Sci. 1989, 267, 941. (6) Doolen, R.; Ray, D. Proc. Soc. Photo.-Opt. Instrum. Eng. 1995, 2547, 364.
Letters (7) Donaldson, D. J.; Guest, J. A.; Goh, M. C. J. Phys. Chem. 1995, 99, 9313. (8) Surface nonlinear spectroscopies have been used by others to study many molecules at liquid interfaces. For a recent review of work in this area, see: Eisenthal, K. B. Chem. ReV. 1996, 96, 1343. (9) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (10) Kirkwood, J. G. J. Chem. Phys. 1930, 3, 300. Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420. (11) Taylor, R. S.; Dang, L. X.; Garrett, B. C. J. Phys Chem. 1996, 100, 11720. (12) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. J. Chem. Phys. 1995, 102, 4574. (13) Matsumoto, M.; Kataoka, Y. J. Chem. Phys. 1988, 88, 3233. Maroncelli, M.; Flemming, G. R. J. Chem. Phys. 1989, 89, 5044. Wilson, M. A.; Pohorille, A.; Pratt, L. R. J. Chem. Phys. 1989, 90, 5211. Townsend, R. M.; Rice, S. A. J. Chem. Phys. 1991, 94, 2207. Benjamin, I. Phys. ReV. Lett. 1994, 73, 2083. (14) Matsumoto, M.; Takaoka, Y.; Kataoka, Y. J. Chem. Phys. 1993, 98, 1464. Benjamin, I.; Pohorille, A. J. Chem. Phys. 1993, 98, 236. Pohorille, A.; Benjamin, I. J. Phys. Chem. 1993, 97, 2664. Chang, T.-M.; Dang, L. X. J. Chem. Phys. 1996, 104, 6772. Tarek, M.; Tobias, D. J.; Klein, M. L. J. Chem. Soc., Faraday Trans. 1996, 92, 559. (15) Huston, S. E.; Rossky, P. J. J. Phys. Chem. 1989, 93, 7888. Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. Dang, L. X. J. Am. Chem. Soc. 1995, 117, 6954. (16) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (17) Taylor, R. S.; Garrett, B. C. The Structure of Alcohols at the Liquid/ Vapor Interface of Water: A Molecular Dynamics Study, in preparation. (18) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016. (19) Pearlman, D. A.; Case, D. A.; Caldwell, J. C.; Seibel, G. L.; Singh, U. C.; Weiner, P.; Kollman, P. A. Amber 4.0; University of California: San Francisco, CA, 1991. (20) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T. A.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; HeadGordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, Revision A.1; Gaussian, Inc.: Pittsburgh, PA, 1995. (21) Besler, B. H.; Merz, Jr., K. M.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431. (22) Berendsen, H. J. C.; Postma, J. P.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (23) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (24) Taylor, R. S.; Schenter, G. K.; Garrett, B. C. Effects of NonEquilibrium SolVation on the Mass Accommodation of Ethanol by H2O; in preparation. (25) Hynes, J. T. In The Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC: Boca Raton, FL, 1985; Vol. 4, p 171. (26) Wilson, J. A.; Pohorille, A. J. Phys. Chem. 1997, 101, 3130. (27) Pratt, L. R.; Pohorille, A. Los Alamos Report Number LA-UR92-3260; Los Alamos National Laboratory: Los Alamos, NM, 1992.
