J. Phys. Chem. B 2000, 104, 4403-4407
4403
Molecular Dynamics Study of Water-Benzene Interactions at the Liquid/Vapor Interface of Water Liem X. Dang* and David Feller Pacific Northwest National Laboratory, P.O. Box 999, Mail Stop K8-91, Richland, Washington 99352 ReceiVed: January 5, 2000; In Final Form: February 15, 2000
The free energy profile of the transport of a benzene molecule across the liquid/vapor interface of water was investigated using molecular dynamics techniques. A large free energy minimum was found near the Gibbs dividing surface. At this free energy minimum, which is an indication of surface activity, the benzene molecule lies parallel to the surface and has an adsorption free energy of about -4 kcal/mol. In addition, there is a free energy barrier to passage from the bulk solution to the surface-adsorbed site. The hydration free energy estimated from the difference between gas-phase and liquid-phase free energies is about -0.1 ( 0.4 kcal/ mol, a value that agrees well with the corresponding experimental measurement of -0.767 kcal/mol. The results described in this paper demonstrate the validity of our approach, as well as the quality of our potential energy surfaces.
I. Introduction The interactions of aromatic hydrocarbon molecules, such as benzene with water, present a special challenge to experimental and theoretical investigations due to the presence of the conjugated π system. These interactions are of great interest as prototypes for the study of hydrophobic hydration interactions, which are ubiquitous in nature. They give rise to the immiscibility of the bulk liquids and play a major role in the formation of micelles and biological membranes.1 Significant experimental and theoretical efforts have focused on studies of benzenewater interactions in solution and in gas-phase clusters. Examples include the extensive experimental work of Zwier and co-workers on the structure and dynamics of benzene-(water)n clusters,2 and quantum chemical calculations on the molecular ice cube-benzene complex by Jordan et al.3 Jonsson, Beveridge, and co-workers4,5 used computer simulation techniques to study the aqueous hydration benzene molecule, and Jorgensen and Severance reported on an extensive study on the hydration of benzene as well as benzene-benzene interactions in water and other solvents using the Monte Carlo method.6 Despite the substantial progress that has been made in this area, a detailed picture of benzene-water interactions has, to our knowledge, not been fully achieved. An accurate treatment of this prototype system is important for developing models capable of handling complex systems in which benzene participates as an environmental contaminant.7 In this paper, we describe the results of classical molecular dynamics (MD) simulations that focus on the mechanisms for transferring a benzene molecule across the liquid/vapor interface of water. Our work is distinguished from earlier contributions by both the computational approach and the potential models. The significance of this work can be summarized as follows: (1) the computed mean force and the corresponding potential of mean force (pmf) allow us to examine the solvent effects on the free energy profile and to infer a mechanism for the transporting processes, (2) adsorption and hydration free energies were estimated from the computed pmf values and can be directly compared with the experimental measurements, and (3) molecular interactions among the species are described using
polarizable potential models, whereas earlier theoretical contributions did not. We show in this paper that the benzene molecule is surface active, i.e., the global minimum in the pmf occurs when the benzene molecule rests at the surface, with the plane of the ring lying parallel to the surface, to maximize the attractive shortrange interaction with water molecules at the interface. The pmf indicates that the benzene is insoluble in water and possesses a very small free energy of hydration (-0.1 kcal/mol), which is in excellent agreement with data obtained from experiment. This work demonstrated the validity of our approaches as well as the water-water and benzene-water potential models. The paper is organized as follows. In section II, we briefly describe the polarizable potential model and computational details used in this study. The simulation results of the benzene molecule moving across the liquid/vapor interface are summarized and discussed in section III. The conclusions are presented in section IV. II. Potential Models and Computational Methods The total intermolecular interaction energy of the system is summarized as follows:
Utot ) Upair + Upol Upair )
∑i ∑j
(1)
( [( ) ( ) ] ) σij
4
12
-
rij
σij
6
+
rij
qiqj
(2)
rij
and N
Upol ) -
µiEi ∑ i)1
0
-
1
N
N
∑∑
2 i)1 j)1,i*j
N
µiTijµj +
|µi|2
∑ i)1 2R
(3) i
Here, rij is the distance between site i and j, q is the charge, and σ and are the Lennard-Jones parameters. E0i is the electric field at site i produced by the fixed charges in the system, µi is the induced dipole moment at atom site i, and Tij is the dipole tensor. The first term in eq 3 represents the
10.1021/jp000054v CCC: $19.00 © 2000 American Chemical Society Published on Web 04/14/2000
4404 J. Phys. Chem. B, Vol. 104, No. 18, 2000
Dang and Feller
attractive charge-dipole interaction. The second term describes the dipole-dipole interaction, and the last term is the energy associated with the generation of the dipole moment µi. During molecular dynamics simulations, a standard iterative selfconsistent field procedure is used to evaluate the induced dipoles. To correctly compute the forces (Coulombic and polarization) using our water model, during the dynamics simulations, the rM, non-atom-centered charge positions were constrained using the following relationship:
rM ) (1 - 2R)rO + R(rH1 + rH2)
(4)
where, rO, rH1 and rH2 are the Cartesian components of the oxygen and hydrogen atom positions. Here R ) rOM/2rOC, rOC is the distance from the oxygen atom to the midpoint of the two hydrogen atoms, and rOM is the distance from the oxygen atom to the M point. The forces FM (Coulombic and polarization) exerted on the M point were calculated at every dynamical time step and redistributed back to the oxygen and hydrogen atoms using the following relationships:
FOW ) (1 - 2R)FM FHW1 ) FH1 + RFM
(5)
FHW2 ) FH2 + RFM FOW, FHW1, and FHW2 are the forces exerted on the oxygen and hydrogen atoms. To avoid unnecessary “polarization catastrophe” due to unfavorable initial conditions, we first equilibrate our system using a nonpolarizable potential model such as TIP4P, and once the system reaches the equilibrium state, we turn on the polarizable potential and continue the equilibration process. The intramolecular potential functions and potential parameters are taken from Amber force fields.8 The rigid-body water-water intermolecular interaction was taken from our earlier work.9 Our water model is rigid, four sites and polarizable, with three point charges similar to the TIP4P format, but we do not consider it to be a modified version of TIP4P. We derived the potential parameters by carrying out extensive molecular dynamics simulations. Our benzene potential model was critically evaluated and found to perform excellently for both the liquid benzene and the equilibrium properties of the liquid/vapor interface of benzene such as the density profile and surface tension. We did not include these results in this paper because we feel that they are irrelevant to this study. The Lennard-Jones parameters for the intermolecular benzene-water potential model were initially constructed using the Lorentz-Berthelot10 combining rules using benzene and water parameters. However, we found that this approach overvalued the binding energy (∼-5.1 kcal/mol) for the waterbenzene interaction and undervalued the distance between the water oxygen and the benzene ring (∼3.0 Å) when compared to accurate, high-level ab initio calculations.11 Therefore, we adjusted the benzene-water Lennard-Jones parameters to correct these deficiencies. These adjustments were based on ab initio potential energy curves which correspond to water approaching the ring from above and from the side. The calculations were performed at the second-order Møller-Plesset perturbation theory (MP2) level with the augmented correlation consistent double-ζ (aug-cc-pVDZ) basis set.12,13 The curves (shown in Figure 1a) were generated by freezing the distance between the oxygen and the benzene center of mass and optimizing all other internal coordinates. The MP2 calculations used the frozen core approximation and were performed with Gaussian 94.14 To
minimize the undesirable effects of basis set superposition error, we corrected the raw potential energy curves with the Boy and Bernardi counterpoise procedure.15 A summary of the benzenewater interaction energy as the water molecule approaches the benzene ring as a function of benzene-water separation is shown in Figure 1a. The final potential parameters used in this work are summarized in Table 1. The benzene atomic polarizabilities were taken from Applequist,16 and the water polarizability was from Stillinger.17 The molecular dynamics simulations are performed on a system consisting of a benzene and 600 water molecules, in a rectangular simulation cell with linear dimensions of 26.22 × 26.22 × 76.22 Å. At the beginning of the simulation, the benzene molecule was located at the center of the rectangular cell. The simulations were performed at a constant volume and temperature with periodic boundary conditions applied in all three directions. The nonbonded interactions (i.e., LennardJones, Coulombic, and polarization) are truncated at molecular center-of mass separations (i.e., water-water and benzenewater interactions were truncated at a molecular separation of 9 Å). During the MD simulations, the SHAKE procedure18 was employed to fix the water geometry by constraining the OH and HH bond lengths. A time step of 2 fs was used in integrating the equations of motion. This approach has been shown to provide results comparable to those of the reaction field or Ewald summation techniques.19 To evaluate the free energies associated with the transfer of a benzene molecule across the vapor/liquid interface, we use the constrained molecular dynamics techniques similar to the approach used for ion transfer across the liquid/liquid interface.20 The reaction coordinate for benzene transfer can be considered as the zs position of the center of mass of the benzene molecule. The Helmholtz free energy difference, ∆F(zs), between a state where the benzene is located at zs, F(zs), and a reference state where the benzene is at z0, F0, is simply
∫zz 〈fz(zs′)〉 dzs′
∆F(zs) ) F(zs) - F0 ) -
s
(4)
0
where fz(zs′) is the z component of the total force exerted on the benzene at a given z position, zs′, averaged over the canonical ensemble. Here, F0 is chosen as the free energy of the system with the benzene located in the bulk water region. During the simulation, the z coordinate of the benzene center of mass was reset to the original value after each dynamical step and the average force acting on the benzene was evaluated. The average forces are subsequently integrated to yield the free energy profile or the potential of mean force (pmf). For the present calculations, the position of the benzene center of mass ranges between z ) 38 Å and z ) 65 Å, with a position increment of 1.0 Å. The total simulation time at each benzene position is at least 350 ps with 50 ps for equilibration. We found the average forces for a given benzene center-of-mass position are converged to within a 100 ps simulation time as demonstrated in Figure 1b. The standard deviation was estimated from the subaverages taken over 10 blocks (30 ps). The average error of the free energy is (0.4 kcal/mol. III. Results and Discussion In Figure 2, we show the pmf for the transfer a benzene molecule across the vapor/liquid interface of water at 298 K. It can be seen from Figure 2 that, as the benzene molecule moves from bulk liquid to the interface, there is existence of a maximum, which indicates that there is a barrier between the bulk solution and the surface-active state. The barrier to
MD Study of Water-Benzene Interactions
J. Phys. Chem. B, Vol. 104, No. 18, 2000 4405
Figure 1. (a, top) MP2(FC)/aug-cc-pVDZ potential curves (with and without counterpoise correction) as a function of the oxygen-benzene center of mass. All other parameters are optimized. The MP2/aug-cc-pVDZ potential curves were shifted so as to match the well depth at the estimated complete basis set limit. The solid circles are the data from the molecular dynamics model. (b, bottom) Cumulative forces (for a given benzene position) acting on the benzene as a function of time.
TABLE 1: Potential Parameters for C6H6, H2O, and C6H6-H2O Interactions Used in the MD Simulationsa molecule
atom type
σ (Å)
(kcal/mol)
q (e)
R (Å3)
C6H6
CB HB HW OW M CB-OW HB-OW CB-HW HB-HW
3.5502 2.5034 0.0000 3.2340 0.000 3.3921 2.8687 2.9377 2.4143
0.0635 0.0250 0.000 0.1825 0.000 0.1077 0.0675 0.0471 0.0296
-0.1430 0.1430 0.5190 0.0000 -1.0380
0.777c 0.172c 0.000 0.000 1.444d
H2Ob
aσ and are the Lennard-Jones parameters, q is the atomic charge, and R is the polarizability. b The potential parameters are taken from ref 9. c Reference 16. d Reference 17.
dissolution from the liquid water is approximately 0.5 kcal/mol. The existence of such a barrier is likely due to the decrease in hydration number of the benzene molecule as it begins to approach the interface. The free energy profile then decreases to a minimum that corresponds to a surface-active state in which the benzene molecule lies on the face of the surface. The location of the minimum is about 1 Å away from the Gibbs dividing surface. The free energy profile then increases as the benzene leaves the interface to become a gas-phase molecule, where the mean force is zero and the corresponding potential of mean force is essentially a constant.
Figure 2. Computed free energy profile of transferring a benzene molecule across the vapor/liquid interface of water.
The hydration free energy estimated from the difference of the gas-phase and the liquid-phase free energies is ∼-0.1 ( 0.4 kcal/mol, in reasonably good agreement with the corresponding experimental measurement of -0.767 kcal/mol.21 The calculated adsorption free energy is -4.0 ( 0.4 kcal/mol. Unfortunately, there are no experimental data that we are aware
4406 J. Phys. Chem. B, Vol. 104, No. 18, 2000
Dang and Feller
Figure 4. The gbo(r) at various benzene center-of-mass positions: 38 Å (solid), 46 Å (dashed), 53 Å (dotted). The inset is a comparison between our simulation and a simulation that used the TIP4P water model.
Figure 3. A snapshot taken at the end of a molecular dynamics simulation when the benzene molecule was located at the liquid/vapor interface of water: (a) side view, (b) top view.
of on the free energy of adsorption of benzene at the water liquid/vapor interface. Eisenthal and co-workers22 investigated properties of phenol at the air/water interface using second harmonic generation techniques and reported a value of -3.8 ( 03 kcal/mol for the adsorption free energy (this is the free energy difference between phenol on the surface and phenol in bulk water). On the basis of the differences between benzene and phenol molecules, we would expect that the experimental value of the benzene adsorption free energy would be larger than -3.8 kcal/mol. To further characterize the benzene transfer free energy across the water vapor/liquid interface, the pmf values were decomposed into energy contributions that correspond to the individual benzene-water interactions. These energy contributions indicated that the computed free energy of benzene transfer is the resultant of a competition between the Lennard-Jones and the electrostatic interactions (Coulombic and polarization). These features are summarized in Figure 2. In Figure 3, a snapshot was taken at the end of a 350 ps segment of the molecular dynamics simulation when the benzene molecule was located at the minimum of the free energy profile. It is clear that the benzene lies parallel to the interface. Radial distribution functions (RDF) and the corresponding coordination numbers provide a more quantitative measure of the local environment of the hydrated benzene molecule in
water. The RDFs for water oxygen atoms relative to the center of mass of the benzene molecule for various benzene positions along the interface (i.e., at the surface-adsorbed site, at the barrier, and in bulk) were computed, and the results are shown in Figure 4. The shoulder in the curve near 3.5 Å represents water molecules situated on the top and bottom of the benzene molecule and coincides with the minimum of the benzenewater potential (see Figure 1a). The peak at 4.9 Å corresponds to the interaction of water molecules along the side of the benzene molecule. Integrating to the first minimum of the gbo(r) yields 30 water molecules in the first hydration shell of the hydrated benzene molecule. At the free energy barrier, the benzene molecule is almost submerged in the liquid as shown by the gbo(r) calculations. The shoulder near 3.5 Å is still visible, and the first hydration shell begins to change. The computed coordination number (n ) 29) is slightly less than the bulk value. The value of the coordination number at the free energy minimum (surface-adsorbed state) is quite small compared to that of the bulk, as one would expect from the gbo(r) results shown in Figure 4. A comparison between our simulation and a simulation that used the TIP4P water model23 as well as the OPLS potential model for benzene also is made in Figure 4. It is clear that both simulations provided nearly identical results with an exception that the first peak position in our gbo(r) calculations shifted a bit to the right. Furthermore, the liquid/ vapor interfaces of water for various benzene positions along the Z axis were examined and found to be similar to that of pure water as the benzene is transported from the aqueous phase to the gas phase. IV. Conclusion In this study, we have established, in detail, a molecularlevel mechanism for transporting benzene across the water vapor/liquid interface, using polarizable potential models and state-of-the art molecular dynamics computer simulation techniques. The results obtained in these studies provided new physical insight into both the free energies and solvent structures as the benzene molecule enters (leaves) liquid water. The benzene molecule is absorbed on the water interface with the plane of the molecule parallel to the water surface. As the benzene leaves the liquid water phase, the first hydration shell of the benzene molecule as a function of the Z axis normal to the interface is significantly reduced. Thus, the mechanism of
MD Study of Water-Benzene Interactions benzene transfer that involves changes in the hydration shell of the benzene is indicated. There is a large free energy minimum near the Gibbs dividing surface, which is an indication that the benzene molecule is surface active. At this minimum, the benzene molecule lies parallel to the surface with an adsorption free energy of about -4 kcal/mol. In addition, there is a free energy barrier to passage from the bulk solution to the surface-adsorbed site. The hydration free energy estimated from the difference between gas-phase and liquid-phase free energies is about -0.1 ( 0.4 kcal/mol, which is quite reasonable agreement with the corresponding experimental measurement of -0.767 kcal/mol. This work demonstrates the validity of our approach as well as the quality of our potential energy surfaces. The differences between the polarizable and effective pair potentials are subjects of great debate, and we have made significant contributions in this area. We have studied the polarization effects on the liquid/vapor interface of water [see the work of Dang and Chang, for example (J. Chem. Phys. 1997, 106, 8149)]. In that study, we have demonstrated clearly that the computed average dipole moments of water molecules near the interface are close to their gas-phase values, while water molecules far from the interface have dipole moments corresponding to their bulk values. Simulations of the vapor/liquid interface that used effective pair potentials, such as TIP4P and SPC/E water models, cannot quantify these effects. It would be interesting to compare our pmf with that obtained from effective pair potentials such as TIP4P and SPC/E water models. We suspect that the pmf obtained using a nonpolarizable potential model would have a characteristic similar to that of our pmf, because our gbo(r) and the gbo(r) using the TIP4P water model are quite similar, as has been demonstrated in this paper. Acknowledgment. This work was performed in the Environmental Molecular Sciences Laboratory at Pacific Northwest National Laboratory under the auspices of the Division of Chemical Sciences, Office of Basic Energy Sciences, U.S. Department of Energy. Battelle operates the Pacific Northwest Laboratory for the Department of Energy. Computer resources were provided by the Division of Chemical Sciences and by the Scientific Computing Staff, Office of Energy Research, at the National Energy Research Supercomputer Center (Berkeley, CA).
J. Phys. Chem. B, Vol. 104, No. 18, 2000 4407 References and Notes (1) Tanford, S. The Hydrophobic Effect; Formation of Miscelles and Biological Membranes, 2nd ed.; Wiley: New York, 1980. (2) Gotch, A.; Zwier, T. S. J. Chem. Phys. 1992, 96, 3388. (3) Gruenloh, C. J.; Carney, J. R.; Arrington, C. A.; Zwier, T. S.; Fredericks, S. Y.; Jordan, K. D. Science 1997, 276, 1678. (4) Linse, P.; Karstrom, G.; Jonnsson, B. J. Am. Chem. Soc. 1984, 106, 4096. (5) Ravishanker, G.; Mehrotra, P. K.; Mezei, M.; Beveridge, D. L. J. Am. Chem. Soc. 1983, 106, 4102. (6) Jorgensen, W. L.; Severance, D. L. J. Am. Chem. Soc. 1990, 112, 4768. (7) Haag, W. R.; Yao, C. C. D. EnViron. Sci. Technol. 1992, 26, 1005. U.S. Department of Energy, Office of Environmental Management, FY 1995 Technology Development Needs Summary, 1994; p 2-17. (8) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Calwell, J. C.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (9) (a) Dang, L. X.; Chang, T.-M. J. Chem. Phys. 1997, 106, 8149. (b) Dang, L. X. J. Phys. Chem. To be submitted. (10) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: London, 1986. (11) Feller, D. J. Phys. Chem. A 1999, 103, 7558. (12) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (13) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. (14) 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.; Andreas, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, E.3; Gaussian, Inc.: Pittsburgh, PA, 1996. (15) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (16) Applequist J. J. Phys. Chem. 1993, 97, 6016. (17) Stillinger, F. H. In The Liquid State of Matter; Fluids, Simple and Complex, Montroll, E. W., Lebowitz, J. L., Eds.; (North-Holland: Amsterdam, 1982; p 341 and references cited therein. (18) Berendsen, H. J. C.; Postma, J. P.; Di Nola, A.; Van Gunsteren, W. F.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. Rykaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (19) Perara, L.; Essmann, U.; Berkowitz, M. J. Chem. Phys. 1995, 102, 450. Sporh, E. J. Chem. Phys. 1997, 107, 6342. (20) Dang, L. X. J. Phys. Chem. B 1999, 103, 8195. (21) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016. (22) Castro, A.; Bhattacharyya; Eisenthal, K. B. J. Chem. Phys. 1991, 95, 8750. (23) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Imprey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926.
