795
J. Phys. Chem. 1986, 90, 795-802 ported oxide. These measurements should be carried out in a zone where NO adsorption does not depend strongly on the temperature (for LaRh03 between 140 and 340 K and P I 50 mmHg) (Figure 8). N O has a number of advantages over other molecules such as C O and O2 used for the characterization of oxides. C O adsorption appears to be more temperature-dependent than NO.29 On the other hand, measurements of adsorption of individual gases and of coadsorption NO-CO seem to indicate that the adsorption bond of C O on as on mixed oxides22is weaker than that of NO. Finally, the method of O2 adsorption at low temperatures (77-195 K) involves the experimental determination of the irreversibly adsorbed part. This is made by the difference of the total and reversible adsorptions. Both these values are, (51) Alexeyev, A.; Terenin, A. J . Cutal. 1965, 4 , 440. (52) Rask6, J.; Solymosi, F. Acta Chim. Acud. Sci. Hung. 1977, 95, 389.
generally, very large and this may constitute an important source of error. Acknowledgment. We are indebted to CAICYT for sponsorship of this work. Thanks are also due to Ministerio de Educacion y Ciencia for a grant awarded to J.M.D.T. and to the SpanishNorth American Joint Committee for Scientific and Technological Cooperation for financial support. Registry No. LaRh03, 11092-15-2; CO, 630-08-0. (53) In this paper the periodic group notation in parentheses is in accord with recent actions by IUPAC and ACS nomenclature committees. A and B notation is eliminated because of wide confusion. Groups IA and IIA become groups 1 and 2. The d-transition elements comprise groups 3 through 12, and the p-block elements comprise groups 13 through 18. (Note that the former Roman number designation is preserved in the last digit of the new numbering: e.g., 111 3 and 13.)
-
Molecular Dynamics Study of the Hydrophobic Interaction in an Aqueous Solution of Krypton Kyoko Watanabet and Hans C. Andersen* Department of Chemistry, Stanford University, Stanford, California 94305 (Received: July 19, 1985)
We have performed molecular dynamics calculations for systems of two and several atomic solutes dissolved in a model for liquid water to study the hydrophobic association of nonpolar solutes in aqueous solution. The parameters of the potentials were chosen to be appropriate for an aqueous solution of Kr. In accordance with previous work, we find that at infinite dilution a pair of atoms has two preferred configurations, one in which they are near neighbors and another in which they are separated by a water molecule, and that the latter is more stable than the former at infinite dilution. In the simulations of systems containing several atomic solutes, the atoms were found to form clusters in which the atomatom distances corresponded to the formation of near-neighbor and solvent-separated pairs. The number of near-neighbor pairs formed is approximately the same (within the noise of the calculation) as there would be in a random distribution of solutes in space at the same concentration, indicating that in water there is no great tendency for hydrophobic association of these solutes. The number of solvent-separated pairs is significantly less than for a random distribution of solutes. On the average, the effect of the solvent is to tend to keep the solutes apart rather than to cause them to associate. This result contradicts the conventional wisdom on hydrophobic interaction, and a new analysis of experimental solubility data gives support to the result.
I. Introduction The phrase “hydrophobic interaction” refers to the supposed tendency for nonpolar solutes and groups dissolved in water to be associated with one another. It is usually assumed to be a major factor influencing the immiscibility of water and nonpolar substances, the thermodynamic properties of aqueous solutions, the conformations of proteins and nucleic acids in solution, association equilibria in aqueous solution, micelle formation by amphiphilic molecules, and the stability of lipid bilayers and biological membranes. The basis for hydrophobic interaction is usually assumed to be the unusual way in which nonpolar molecules and groups are hydrated in water. This unusual hydration is often called “hydrophobic hydration”. Computer simulation studies have contributed greatly to our understanding of the hydrophobic interaction. Pangali et a1.I calculated the potential of mean force, or equivalently the solute-solute correlation function, for a pair of nonpolar solutes in water using their force-bias Monte Carlo technique. Their results indicated the existence of two relatively stable configurations for the solute atoms, in agreement with the earlier semiempirical theory of Pratt and Chandler.2 One of these configurations corresponds to two solutes in contact. The interesting feature of the results was the greater stability of the second configuration, in which two solute atoms are separated by a water molecule. The Natural Sciences and Engineering Research Council of Canada postdoctoral fellow.
0022-3654 I86 12090-0795$01.5010
existence of two stable configurations was also implied by two other computer simulation s t u d i e ~ ,which ~ . ~ calculated the mean relative force between nonpolar solutes in water. Further insights into the problem of hydrophobic interaction were provided by molecular dynamics studies in which the motion of solutes in water was calculated. Geiger et aL5 studied a system consisting of 2 Lennard-Jones solute atoms and 214 water molecules. The solute atoms, which were initially in contact, drifted to permit a single water molecule between them, resulting in a formation of vertex-sharing solvent cages containing the solutes. Using a similar approach, Rapoport and Scheraga6 calculated the motions in a system composed of 4 solute atoms and 339 water molecules. The simulation was run for 70 ps, which was considerably longer than the run time of 8 ps used by Geiger et al. This simulation, however, failed to reveal any tendency for the initially well-separated solutes to associate. As pointed out by Rapoport and Scheraga, however, there is a possibility that the duration of their simulation was insufficient to permit observation of the association of solutes. (1) C. Pangali, M. Rao, and B. J. Berne, J . Chem. Phys., 71, 2975 (1979). (2) L. R. Pratt and D. Chandler, J . Chem. Phys., 67, 3683 (1977). (3) C. S. Pangali, M. Rao, and B. J. Berne in “Computer Modeling of Matter“, P. Lykos, Ed., American Chemical Society, Washington, D.C., 1978. (4) S.Swaminathan and D. L. Beveridge, J . A m . Chem. Soc., 101, 5832 (1979). ( 5 ) A. Geiger, A. Rahman, and F. H. Stillinger, J . Chem. Phys., 70, 263 ( 1979). (6) D. C. R a p p o r t and H. A. Scheraga, J . Phys. Chem., 86, 873 (1982).
0 1986 American Chemical Societv
796 The Journal of Physical Chemistry, Vol. 90, No. 5, 1986
There have also been several computer simulation studies of single hydrophobic solutes in aqueous solution to learn about the thermodynamic and structural properties of such solutions at infinite The solution theory of McMillan and Mayer2z provides a convenient way of thinking about the distribution of solutes in a solution. According to this theory, the solutes in a solution are distributed in the same way as in a hypothetical gas containing only the solutes. This gas has the same concentration of solutes as does the solution, and the interaction potentials in the hypothetical gas are related to the solute potentials of mean force at infinite dilution in the solution. The potential in the hypothetical gas is N
N
Watanabe and Andersen
-2t I
Each of the wn(0)becomes zero in the limit that one or more of the n particles become separated from the others. The w,(O) are the potentials of mean force at infinite dilution. They can in principle be calculated from the distributions of 2, 3, ..., n solutes in solution at infinite dilution. Once they are known, they can in principle be used to calculate the distribution of solutes at finite concentration. A significant question for understanding a solution is that of whether the potentials are pairwise additive, i.e., whether only the w2(O)terms need to be kept for an accurate description of the solution. This paper reports the results of molecular dynamics computer simulation studies of hydrophobic interaction in an aqueous solution of Kr both at infinite dilution and at finite concentrations of about 0.8 and 1.4 M. The results at infinite dilution confirm those of Pangali et al.' The solute-solute pair correlation function at finite concentrations is found to be consistent with the pairwise additivity of the potential of mean force at infinite dilution. The systems simulated exhibited no strong tendency for the solutes to associate. When compared with randomly distributed particles, the nonpolar atoms dissolved in water were found to be correlated in such a way that they stayed away from each other. These results are in conflict with the usual notions of hydrophobic association of nonpolar solutes. A new analysis of published data on the solubility of nonpolar gases in water gives support to these results. 11. Intermolecular Forces and Simulation Method The water-water intermolecular potential used in this work is a slight modification of the revised central force model of Stillinger and RahmanZ3 The modifications were a truncation of the potential at large distances using a switching function and a slight modification of the oxygen-oxygen potential so that the model gave approximately the correct density at a pressure of 1 atm at 29.5 O C . The potential is discussed in more detail by Swope and A n d e r ~ e n . ~The ~ water intramolecular potential was a set of (7) G. Dashevsky and G. N. Sarkisov, Mol. Phys., 27, 1271 (1974). (8) J. C. Owicki and H. A. Scheraga, J . Am. Chem. Soc., 99, 7413 (1977). (9) S. Swaminathan, S . W. Harrison, and D. L. Beveridge, J . Am. Chem. SOC.,100, 5705 (1978). (10) S . Okazaki, K. Nakanishi, H. Touhara, and Y . Adachi, J . Chem. Phys., 71, 2421 (1979). (11) P. J. Rossky and M. Karplus, J . Am. Chem. Soc., 101, 1913 (1979). (12) G. Alagona and A. Tani, J . Chem. Phys., 72, 580 (1980). (13) S . Okazaki, K. Nakanishi, H. Touhara, and Y . Adachi, J . Chem. Phys., 72, 4253 (1980). (14) G. Bolis and E. Clementi, Chem. Phys. Lett., 82, 147 (1981). (15) K. Nakanishi, S . Okazaki, K. Ikari, and H. Touhara, Chem. Phys. Lett.. 84. 428 (1981). (16) S . OkizakilJ. Nakanishi, H. Touhara, and N. Watanabe, J . Chem. Phys., 74, 5863 (1981). (17) G. Bolis, G. Corongiu, and E. Clementi, Chem. Phys. Lett., 86, 299 119821 \----,. (18) G. Alagona and A. Tani, Chem. Phys. Lett., 87, 337 (1982). (19) W. L. Jorgensen, J . Chem. Phys., 77, 5757 (1982). (20) R. Rosenberg, R. Mikkilineni, and B. J. Berne, J . Am. Chem. SOC., 104, 7647 (1982). (21) W . L. Jorgensen and J. D. Madura, J . Am. Chem. SOC.,105, 1407 (1983). (22) W. G. McMillan and J. E. Mayer, J . Chem. Phys., 13, 276 (1945). (23) F. H. Stillinger and A. Rahman, J . Chem. Phys., 68, 666 (1978). ~~
1
6
6
10 r
12
14
I
(A)
Figure 1. Mean force, ( F ( r ) ) ,between two solutes dissolved in water. The points are molecular dynamics results for infinite dilution: circles, 64 water molecules; triangles, 148 water molecules; squares, 198 water molecules. The smooth curve was obtained by analytic differentiation of a function that was fit to the finite-concentration results (6 solutes in 394 water molecules) for the pair correlation function.
simple atom-atom harmonic potentials as used by Andrea et al.25 The watersolute potential was assumed to be a Lennard- Jones interaction between the solute atom and the oxygen of the water molecule. Swope and A n d e r ~ e nhave ~ ~ calculated the solubility of such solutes in this model of water as a function of temperature for a range of values of the Lennard-Jones parameters for the atom-water interaction. The values chosen in this work, e = 0.395 kcal/mol and u = 3.44 A, led to a solubility with about the same magnitude and temperature dependence as that observed experimentally for Kr. The solutesolute potential was assumed to be a Lennard-Jones interaction with E = 0.401 kcal/mol and u = 3.6 A. These parameters were chosen to make the potential similar to the Kr-Kr potential of Barker et aLZ6 The molecular dynamics calculations were performed for a cubic system of constant volume using periodic boundary conditions. The velocity form of the Verlet algorithmz7 was used to integrate Newton's equations of motion with a time step of 0.6 fs. Unless otherwise noted, the calculations used the constant-temperature molecular dynamics method;**stochastic collisions corresponding to a bath temperature of 298 K were applied to all the atoms and molecules every 0.3 or 0.6 ps. 111. Hydrophobic Interaction of Two Solutes For the case of two solute atoms dissolved in water, the potential of mean force, wz(r), and the solute-solute pair correlation function, gss(r), are related by
where k is the Boltzmann constant and T i s the absolute temperature. (Note that the w2 in this equation is the potential of mean force at the concentration of interest, not the potential of mean force at infinite dilution, which is denoted w2(0).) If two solutes are fixed at a distance r from one another, the force on each atom, averaged over a canonical distribution for the solvent, denoted (F(r)), is related to w(r) in the following way (F(r)) = -dwz(r)/dr
(2)
The direction of (F(r)) for an atom is along the intersolute vector, with positive values of ( F ( r ) )denoting that a solute is on the (24) W. C. Swope and H. C. Andersen, J . Phys. Chem., 88,6548 (1984). (25) T. A. Andrea, W. C. Swope, and H. C. Andersen, J . Chem. Phys., 79, 4576 (1983). (26) J. A. Barker, R. 0. Watts, J. K. Lee, T. P. Schafer, and Y . T. Lee, J . Chem. Phys., 61, 3081 (1974). (27) W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, J . Chem. Phys., 76, 637 (1982). (28) H. C. Andersen, J . Chem. Phys., 72, 2384 (1980).
Hydrophobic Interaction in an Aqueous Kr Solution average being repelled from the other solute. ( F ( r ) ) can be calculated by computer simulations by keeping the two solutes fixed in position while allowing the solvent molecules to move freely. This method, used by Pangali et a L 3 is the one we used in this work. Simulations were performed for systems containing two solute atoms and various numbers of water molecules to evaluate ( F ( r ) ) for a range of intersolute distances. Systems containing 64, 148, and 198 water molecules were used for calculations at intersolute distances ranging from 3.4 to 6.9 A, 5.8 to 11.2 A, and 9.6 to 14.4 A, respectively. For each simulation the system volume was determined from the molar volume of Kr at its boiling point, 34.7 ~ m ~ / m oand l , ~the ~ pure-water density appropriate to the water model used in this work. The duration of each run was between 180 and 450 ps. The distances between the solute atoms and their images in the neighboring boxes were maximized by placing the solutes along the diagonal direction of the box. The average force on each atom in the direction of the intersolute axis as well as the average force in the two perpendicular directions was computed for each value of time and averaged over time to obtain the desired ensemble averages. Because of the symmetry of the system, the ensemble average components of the forces in directions perpendicular to the intersolute coordinate should be zero even for a system of finite size. The results are shown in Figure 1. The ( F ( r ) ) given there is actually the average of the results for both solutes. The statistical error was estimated from the magnitudes of the four time-averaged components of the force perpendicular to the intersolute axis. The oscillatory nature of the solvent-averaged force as a function of the intersolute distance is evident. The intersolute distances at which the force vanishes correspond to the extrema in the wz(r) and gss(r),with the smallest such distance occurring close to the u value for the solutesolute interaction potential. This distance corresponds to that for a solute pair in a state of near-neighbor hydrophobic association. The second maximum in gsj(r)falls at r = 7 A, which is nearly twice the u value for the solute-water interaction, and corresponds to the solvent-separated hydrophobic association of a pair of solutes. The third maximum in gss(r)may be located at 10 A, a solute-solute distance probably corresponding to two solute atoms separated by two water molecules. Presence of extrema at distances greater than r 11 A cannot be confirmed owing to the magnitudes of the statistical errors associated with evaluation of ( F ( r ) ) . Use of the calculated forces in eq 1 and 2 in principle enables the function gjj(r)to be calculated. However, the magnitudes of the errors in these forces lead to large uncertainties in the derived function. There is an additional source of uncertainty due to an undetermined constant of integration in wZ(r),which leads to an unknown multiplicative constant in gSs(r).For these reasons, we do not present the wz(r)and g,(r) determined using the simulation results for ( F ( r ) ) . These infinite-dilution results will be compared below with results for finite concentrations. Although the construction of gss(r)from the mean force data is subject to these uncertainties, various ways of doing this with our data always led to the conclusion that solvent-separated pairs are favored over near-neighbor pairs, as found in previous theories.'q2
-
-
IV. Hydrophobic Association of Several Solutes Hydrophobic association of several nonpolar solutes was studied by simulating a system of 5 solute atoms and 195 water molecules. The initial configuration was generated from an equilibrated system of 200 water molecules, in a volume corresponding to a density of 1 .O g/cm3, by randomly selecting five water molecules and replacing them by solute atoms located at the former positions of the oxygen atoms. Equilibration was carried out for 30 ps with long-ranged purely repulsive interactions acting between the solutes while the water-water and solute-water intermolecular forces remain as described in section 11. This part of the simulation (29) K. W. Miller and J. H. Hildebrand, J . Am. Chem. SOC.,90, 3001 (1968).
The Journal of Physical Chemistry, Vol. 90, No. 5, 1986 797 12 2
0 O
@ 0
168
120
156
192
264
1
276
384
420
I
456
492
516
Q
@-a 0 0
/ b
I
0 Q
600
v
O
0
,@
816
1152
1020 0
912
9 1440
0
1236
0 @
0 1536
i Figure 2. Solute configuration for five solutes initially well separated in a system containing 195 water molecules at (a) t = 0-480 ps and (b) t > 480 ps. A solid line between two particles denotes that they are a near-neighbor pair, with interatomic distance less than 5.2 A. A dashed line denotes that they are a solvent-separated pair, with interatomic distance between 5.2 and 8.6 A. The number by each diagram is the time in picoseconds.
permits the solute atoms to separate from each other and allows the water molecules to reorganize around each solute atom. At the end of the equilibration time, the solute-solute intermolecular force was changed to the Lennard-Jones form given in section 11. This time is chosen as t = 0. The simulation was continued, and the state of the system was observed at 12-ps intervals. The solute configuration is shown schematically in Figure 2a,b. In describing the configuration of a pair of solute atoms, we will use the following conventions: (i) Two solutes will be referred to as a near-neighbor pair if the intersolute distance is less than 5.2 A, the separation distance corresponding to the first minimum in the solute-solute pair correlation function at infinite dilution.
Watanabe and Andersen
798 The Journal of Physical Chemistry, Vol. 90, No. 5, 1986 (ii) Solute pairs with separation distances in the 5.2-8.6-A range will be referred to as solvent-separated pairs. (iii) Solute pairs separated by more than 8.6 A, an estimated intersolute distance for the second minimum in the solute-solute pair correlation function at infinite dilution, will be considered unassociated. At t = 0, the minimum image interatomic distances for the five solutes ranged from 8.9 to 12.9 A. There were no near-neighbor pairs or solvent-separated pairs. At t = 12 ps, a cluster of four particles containing no near-neighbor pairs had formed. The cluster spanned the periodic cell of the simulation. At t = I20 ps, the cluster had grown to include all five particles. The fiveparticle cluster was compact, and it did not span the periodic cell. This compact five-particle cluster persisted to t = 264 ps. From t = 276 to 384 ps, the five-particle cluster spanned the periodic cell, except at t = 324 ps. The association of all five particles continued for an additional 100 ps to t = 480 ps. Typical configurations observed during this period are shown in Figure 2a. The clusters formed during the first 480 ps of the simulation contained primarily solvent-separated pairs rather than nearneighbor pairs of solutes. The five-particle cluster started to break down at t = 492 ps. A significant decrease in the degree of association among the solute atoms was observed around this time. The simulation was continued for approximately one more nanosecond. During this period there was no tendency for the development of the relatively strong solute associations that were observed during the initial 0.5 ns of the simulation. During the last nanosecond of the simulation, the solute atoms formed compact clusters involving two to five atoms, with clusters extending across the periodic cell appearing only transiently. Typical solute configurations observed during this period of the simulation are shown in Figure 2b. One measure of the extent of aggregation of the solutes is the average number of near-neighbor pairs formed. Let nnn denote this quantity. Similarly, let ns denote the average number of solute pairs at solvent-separated distances. For the time period from t = 468 to t = 1536 ps, nnnand nsswere 0.90 0.14 and 2.85 f 0.26, respectively. (The statistical errors were estimated by breaking the long simulation into 120-ps intervals and assuming that each 120-ps interval is statistically independent of the others.) For comparison, if the five solutes were randomly distributed in the box, the numbers would be 0.679 and 3.441, respectively. The overall hydrophobic association for this period is significantly smaller than that observed from t = 12 to t = 456 ps, during which the average numbers of pairs were 0.53 f 0.20 and 5.63 i 0.17 for the near-neighbor and solvent-separated association, respectively. Because of the qualitative and quantitative differences between the results of the first half-nanosecond and those of the last nanosecond, we interpret the configurations in the last nanosecond as being the true equilibrium states and the cluster seen in the first half-nanosecond as a long-lived transient. The results of this simulation give no support to the conventional wisdom of hydrophobic interaction that suggests that the solutes should want to aggregate and form several near-neighbor pairs. It is still possible, however, that given a sufficiently long time they might form such a stable cluster and that the simulation was not long enough for this to happen. To test whether a compact cluster with many near-neighbor pairs is actually favored in water, we performed a second simulation in which the solutes initially formed such a compact cluster and we observed whether the structure was stable and persisted. The initial configuration for this simulation was obtained by constructing a five-atom cluster containing several near-neighbor pairs. The water molecules were equilibrated for 23 ps with the solute atoms held fixed. The velocities were rescaled periodically during this equilibation to maintain a temperature of 298 K. Further calculations were carried out using the same procedure as in the earlier simulation involving a system of 5 solute atoms and 195 water molecules. Solute configurations were observed at 12-ps intervals as before. The solute configuration is shown schematically in Figure 3. The five solutes remained closely associated for approximately 180 ps with some configurational rearrangements. The cluster remained compact with five to eight near-neighbor pairs, and every
*
120
I
228p240
I
60
@-q
3 4
180
192
8@& 2 -.....
,A
- 4
252
264
288
312
336
348
8
0
o b
0 Figure 3. Solute configuration for five solutes initially forming a compact cluster in a system containing 195 water molecules. (See caption of Figure 2 for the meaning of the figures.)
solute formed a near-neighbor pair or a solvent-separated pair with every other solute most of the time. The configuration at t = 180 ps contained six near-neighbor and three solvent-separated pairs. Then the cluster became less closely associated. For t > 180 ps, there are at most four near-neighbor pairs, and for t > 252 ps there are at most two. In the last part of the simulation, from t = 228 to I = 348 ps, the structure was similar to that observed in the last nanosecond of the simulation in which the solute atoms were well separated initially. We conclude that the initial compact cluster is unstable in water. Thus, starting from both an initial condition in which the atoms in water are well separated and an initial condition in which they form a cluster with several near-neighbor pairs, we found the system evolved toward states in which the atoms are not highly associated and in which the association is primarily through solvent-separated pairs. Since the same sorts of final configurations were achieved in both cases, we conclude that the runs were long enough to give a true indication of the equilibrium states. The systems in the two previous simulations were small. They contained 195 water molecules, and when the five atoms formed a cluster containing primarily solvent-separated intersolute distances, the cluster could span the periodic simulation cell. The small system size and periodic boundary conditions might be creating artifacts. To check this, we performed a third simulation for a larger system containing 6 solute atoms and 394 water molecules (nominal solute concentration of about 0.8 M). The volume of the new system was twice that used for the five-solute simulations. After thorough equilibration, the solute configuration was observed every 12 ps as before. In addition, gss(r)was calculated from solute coordinates taken every 0.06 ps. An observation of the solute configuration as a function of time indicated that the effect of the initial configuration persisted for a long time as before. The values of nnnand nss for the initial 468 ps were 0.45 f 0.1 1 and 2.98 0.1 7 , respectively. Corresponding numbers obtained from t = 480 ps to the end of the simulation at t = 1068 ps were 0.66 f 0.19 and 2.02 f 0.12, respectively. Figure 4 shows some of the solute configurations observed during the last 0.6 ns of the simulation. As in the previous simulation, we regard the initial half-picosecond as an initial transient, and the last half of the simulation is regarded as true equilibrium states. Figure 5 shows the gss(r)calculated from the last 0.6 ns of the simulation of 6 atoms in 394 water molecules. Also shown in the figure is the g,(r) obtained in the last nanosecond of the simulation
*
The Journal of Physical Chemistry, Vol. 90, No. 5. 1986 799
Hydrophobic Interaction in an Aqueous Kr Solution
3
576
4
@ o
0
0 1 I
624
1
9,s
720
672 0 0
1
@
@
c
864
0
0
4
6
8
r
- Q
10
(A)
Figure 5. Solutesolute pair correlation function, g,(r), determined from
simulation of system containing several solutes in water. The circles represent values calculated for the system containing 6 solutes and 394 water molecules. The filled circles represent values determined for the system containing 5 solute atoms and 195 water molecules.
Figure 4. Solute configuration for six solutes in a system containing 394 water molecules. (See caption of Figure 2 for the meaning of the figures.) TABLE I: Average Number of Near-Neighbor and Solvent-Separated Pairs
simulation
nnn
n,,
random distribution of 6 atomsb
0.90 & 0.14 0.679 0.66 f 0.19 0.509
2.85 f 0.26 3.441 2.02 f 0.12 2.580
+ 195H,0a random distribution of 5 atomsb 6Kr + 394HzOc 5Kr
Calculated for the last nanosecond of the 1.5-11s molecular dynamics run. *Analytic calculation. CCalculatedfor the last 0.6 ns of the 1.1-ns molecular dynamics run. of the system containing 5 solute atoms and 195 water molecules. The two results for gS(r) are in good agreement despite the fact that they were obtained for two different concentrations. Table I summarizes the results for nnn and nss for the later equilibrated parts of first and third simulations. For comparison, the table gives the numbers that correspond to a completely random distribution of the solutes in the volume. For both concentrations, nnn, the average number of near-neighbor pairs, is slightly greater than for a random distribution, but the excess is of the order of the statistical noise of the calculation (Le., the excess is about 1.6 and 0.8, respectively, of the standard errors for calculated values of nnnat the two concentrations). In Figure 5 this corresponds to the fact that in the vicinity of the first peak of g,, the curve is greater than one for some distances and less than one for others, and the net number of near-neighbor pairs is about the same as for a random distribution. On the other hand, for both concentrations the average number of solvent-separated pairs is less than that for a random distribution and the difference appears to be significant (Le., the difference is about 2.3 and 4.7 times the standard error of n,,). In Figure 5 this corresponds to the fact that g,, in the vicinity of the second-neighbor peak is consistently less than unity. The sum of nnnand nss for both concentrations is significantly less than for a random distribution. As noted above, runs at the two concentrations led to the same results for gs(r). This suggests that both concentrations are low enough that the following low-concentration limiting results hold: w2(r) = w2(0)(r) g,(r) = exP(-w2'0'(r)/kT)
(3b)
(The right side of the second equation is independent of con-
centration.) Two additional checks were performed to verify this. First, this relationship suggests that the finite concentrations g,, can be calculated from the infinite-dilution potential of mean force. In section I11 we presented results for the mean force at infinite dilution. If those results could be integrated over distance, the potential of mean force at infinite dilution could be constructed and g, calculated. In practice this is difficult because of the noise in the mean force calculations and the unknown constant of integration. As an alternative, we can work in the opposite direction. The finite-concentration g,,(r) results were fit by a continuous differentiable function. The logarithm of this yields the potential of mean force at finite concentration, and the derivative of the latter gives the mean force at finite concentration. The results are plotted in Figure 1 as the smooth curve. The agreement of the curve with the calculated infinite-dilution points confirms the validity of eq 3 at these concentrations. An additional check was performed by constructing hypothetical potentials of mean force at infinite dilution that were consistent with the infinite-dilution mean force results and that contained arbitrary choices of the unknown constant of integration. Simulations of gases of atoms with those potentials at the atomic concentration of the solution simulations led to pair correlation functions consistent with eq 3. A corollary of this discussion is that although the potentials of mean force at infinite dilution may not be pairwise additive (Le., the w3(0) and subsequent terms in eq 1 may be nonzero), the nonadditivity has no effect on g,, at these concentrations. A striking consequence of Figure 5 is that on the average the solutes in solution tend to stay away from each other. They exhibit what might facetiously be called "hydrophobic repulsion" rather than hydrophobic attraction. This suggests that although they do not want to go into the water (Le., the solubility is low), and hence they are in a sense hydrophobic, once they are in solution they prefer to be surrounded by water, and in that sense they are hydrophilic. This tendency to avoid one another in solution is described quantitatively by the osmotic second virial coefficient, given by B, = ( 1 / 2 ) x m [ 1 - g,,(r)]4ar2 dr In evaluating the integral, we assumed the g,,(r) to be unity for r greater than 11.4 A. We obtained a large positive value, 502 A). The correlation function to be used in this equation should be the solute-solute correlation function at infinite dilution. We used our finite-concentration results, but the comparisons above demonstrate that the results are independent of concentration and are hence probably close to the infinite-dilution results. For future reference, this also suggests that the nonideality of aqueous solutions of Kr is dominated by the second osmotic virial coefficient contribution at these concentrations. The Appendix contains an analysis of data on the pressure dependence of several small nonpolar gaseous solutes in water.
800 The Journal of Physical Chemistry, Vol. 90, No. 5, 1986 12 0
Watanabe and Andersen TABLE II: Osmotic Second Virial Coefficients Calculated from Solubility Measurements of Gases in Water
temp. He
Ar
37.8 0.2
H2
B,.
A'
-784
CH4
NZ
.5
O C
25 25 50 0 25 50 25
-525 -335 -164 -300 -194 -33 -1 18 459
25 "C
0
0.002
0004
0006
000,
xs
Figure 6. Analysis of solubility data of gases in water to obtain Henry's law constants (from the intercepts) and osmotic second virial coefficients from the slopes. See eq S,6, and 9 in the Appendix. Straight lines are the results of thi: least-squares fit to the experimental data after eliminating some data as explained in the text.
That analysis confirms that at these concentrations the second osmotic coefficient dominates the solution nonideality and suggests that it is likely that Kr does have a large positive value of B,.
V. Summary and Discussion In this work we have used intermolecular potentials that are consistent with some of the known properties of water and aqueous solutions of krypton. The water-water potential is a slight modification of the revised central force model, which was constructed to fit many of the properties of pure water. The solute-water interaction was one that has been shown to give temperature-dependent solubilities in agreement with experiment for Kr. The solute-solute potential is similar to an accurate Kr-Kr potential. It is important to use potentials that have some degree of quantitative validity as well as qualitative validity when investigating hydrophobic effects, because the work of Pratt and Chandler2s30indicates that such effects are sensitive to the relative strengths of the various interactions in the solution. We find that krypton atoms in water tend to stay away from each other rather than to aggregate; Le., they have a large positive value of the second osmotic virial coefficient. When they are close together, they tend to form near-neighbor or solvent-separated pairs, with the latter being more favored than the former, as has been found by previous computer simulations, but the total number of such close encounters is less than for a random distribution of s o h tes. These results are consistent with a new analysis of solubility data for small nonpolar solutes. They are also consistent with the theorv of Pratt and Chandler.30 This theory predicts that if solute-solute attractions are much more important than solute-solvent attractions, then solutes will tend to cluster together in water. If, on the other hand, solutesolvent attractions are more important, the solutes will tend to be surrounded by water and will not want to cluster together. If it is hypothesized that the relative strer gths of the solventsolute and solutesolute attractions change in tlie homologous series He, Ne, Ar, Kr, Xe in such a way that solute-solvent attractions become more important at the end of the series, then the larger solutes should indeed stay away from each other in solution. In fact, the larger solubility of the (30) L. R. Pratt and D. Chandler, J . Chem. Phys., 73, 3434 (1980)
In K,
Figure 7. A correlation between the Henry's law constant and the osmotic second virial coefficient for small nonpolar solutes in water at 25 OC. The circles indicate values determined from the solubility measurements. The Henry's law constants for Ne, Ar, Kr, and Xe used for the extrapolation are determined from the solubility data at 25 'C and 1 atm listed by Miller and H i l d e b r a ~ ~ d . ~ ~
larger solutes is probably due to the larger strength of their attractions to water, which is consistent with this hypothesis. Our results clearly point out the necessity for performing long runs if phenomena of this type are to be investigated by computer simulation. To obtain accurate mean forces, the runs must be many multiples of 30 ps long, since successive 30-ps runs gave significantly different results. To investigate hydrophobic association, long runs are needed. The compact cluster of five atoms with several near-neighbor bonds had a lifetime of about 180 ps, but it does not represent the dominant equilibrium configuration for the atoms. The more diffuse cluster found in the initial simulation fluctuated in shape and had a lifetime of about 500 ps, but it too does not represent the dominant equilibrium configuration. In these simulations, runs of too short a duration would give quantitative results that are subject to large statistical error and qualitative results about the absolute stability of configurations that are incorrect.
Acknowledgment. This work was supported by the National Science Foundation (Grants CHE81-07165 and CHE84-10701). Appendix The osmotic second virial coefficient may in principle be determined from measurements of deviations from the thermodynamics of the ideal dilute solution. One example of the ideal behavior of solutions is Henry's law for the solubilities of gases in liquids in the limit of infinite dilution, according to which there exists a direct proportionality between the mole fraction of solute in solution and its vapor pressure in the vapor phase. In the following an expression will be derived for deviations from Henry's law at finite solute concentrations, and available solubility measurements of nonpolar molecules in water together with measurements of nonideal behavior of the gases will be used to calculate osmotic second virial coefficients of the solutions. The system of interest contains two components, solute and solvent molecules, and consists of a liquid phase in equilibrium with a gas phase. We consider the special case where the solubility of solute is small and the vapor pressure of solvent is small compared to the total pressure. Clearly, the solubilities of nonpolar
The Journal of Physical Chemistry, Vol. 90, No. 5, 1986 801
Hydrophobic Interaction in an Aqueous Kr Solution molecules in water are of this type. When the vapor pressure of solvent is small, the total pressure, P, can be identified with the partial pressure of the solute and the vapor phase can be treated as consisting of solute molecules alone. The chemical potential, pg, of the solute molecules in the vapor phase is then given by
kL,(T,P)= P ( T ) + k T In I f ( T , P ) / P l where p e ( T ) is the standard chemical potential of the solute molecules in the gas phase at temperature T and f ( T , P ) is the fugacity. Here P denotes the standard pressure of 1 atm. The chemical potential, p,, of the solute molecules in the liquid phase can be written as
r,(T,P,x,) = k T In x, + D,(T)
+ D2(T)(P- P ) + D 3 ( T ) x ,+ O( ),
(4)
where x , is the mole fraction of solute molecules in the solution and O( ), represents terms of second order or higher in (P - P ) or x,. D2(T ) and D3( T), which are thermodynamic quantities related to derivatives of the chemical potential at infinite dilution, can be evaluated by standard thermodynamic reasoning to give DdT) = vm,O
kT D3(T) = ~ ( 2 B +2',u L."w
- 2u,:
(5)
+ kTK,')
(6)
where u,: is the partial molar volume of solute in solution, umWo is the partial molar volume of solvent, and K', is the isothermal compressibility, evaluated at infinite dilution. The equality of the chemical potential in the liquid and vapor phases gives
At low pressures this reduces to c
If one identifies the right side of eq 8 as the logarithm of the Henry's law constant, KH,eq 7 becomes
In
[
= In KH
1
+ E[D,(T)P + D,(T)x,] + O( )* (9)
The only approximations made in deriving this result are that the vapor pressure of the solvent is negligible compared to the pressure and that the solution is dilute enough that quadratic and higher order terms in eq 4 for the chemical potential of the solute are negligible. In particular, no assumption is made that the gas phase of the solute is ideal. This result should be applicable then in situations in which the deviations from ideal behavior of the solute are dominated by the second osmotic virial coefficient. Thus, solubility data at high pressures can be analyzed by plotting In ( f / x , P ) - umSeP/kT,which will be denoted by f in the remainder of the discussion, against x,. T h e y intercept of such a plot gives In KH,and the low-pressure limiting slope, which is equal to D 3 / k T , enables B2 to be determined from eq 6. We have analyzed the solubility measurements in water at high pressures for He,31H,,32 N2,3,CH4,33and Ar.34 It should be noted that the accuracy of B2 determined from the present method depends on the accuracy of solubility measurements as well as (31) R. Wiebe and V. L. Gaddy, J . A m . Chem. Sor., 57, 847 (1935). (32) I. R. Krichevsky and J. S . Kasamovsky, J . Am. Chem. SOC.,57, 2168 (1935). (33) 0. L. Culberson, Trans. A m . Inst. Min., Metaii. Pet. Eng., 192, 223 (1951). (34) B. Sisskind and I. Kasarnowsky, Z . Anorg. Chem., 200, 279 (1931).
on the accuracy of the values used for.:,v We used the experimentally determined v a l ~ e s ~of ~',-v ~ ' for all of the substances considered here except for He. Due to the apparent lack of experimental measurement on: ,0 for He, an approximation was . ~ Ar, ~ made using the van der Waals b value, 23.7 ~ m ~ / m o lFor we used the experimentally determined value appropriate to 25 "C and applied it to the analysis of the solubility data for 0.2 "C. We expect the effects of these approximations on the calculated B2 to be not serious. The values of the fugacities of N2 and H, are those of Deming and S h ~ p e . ~ ~For . ~ other ' substances, the fugacities were calculated from
by using the compressibility factor Z given in ref 41, 42, and 43 for He, CH4 and Ar, respectively. The values o f f calculated from these data are plotted against x, in Figure 6. It can be seen that in general the solubility measurements show a good internal consistency. Exceptions are observed for the lowest concentration data in the isotherm of CHI at 25 "C and that of Ar. We note that the concentrations at which these inconsistenciesoccur ( x , 0.0005 and 0.001) are well inside the concentration range where the linearity between f and x , is exhibited by the other isotherms shown in Figure 6. Furthermore, for the two isotherms in question it can be seen that linearity holds very well for the higher concentration data points. For these reasons it seems reasonable to assume these inconsistencies to be due to experimental errors and to ignore these two points in determining the limiting slopes of the isotherms. The straight lines drawn in Figure 6, except for the N2 isotherms, were obtained from the least-squares fit to the data shown in the figure, after elimination of the two data points discussed above. For isotherms of N2 solutions, the deviations from the linearity seem to start at around x, = 0.003, and therefore straight lines were fitted to data points in the range x , < 0.003. The values for B2 calculated from these straight lines are summarized in Table 11. One observation that can be made from Figure 6 is the increase in the slope of the isotherms as their ordinate values decrease, indicating an increase in the values of B2 as the solubilities of gases in water increase. This trend is confirmed in Figure 7 where B2 is plotted against In KH for He, N2, H2, and CH,. A qualitative prediction of the values of B2 for other substances was attempted by fitting a straight line to the experimental values for He through CH4 (see Figure 7). The values of B (at 25 "C) predicted in this for Ne, Ar, Kr, and Xe, way are -682, -13, 301, and 644 respectively. The values of B2 given in Table I1 exhibit a significant temperature dependence. In particular, the values calculated for H, indicate that B2 attains its minimum value at approximately 25 "C. For CH4 the minimum in B, appears to occur at a higher temperature. In view of such temperature dependence of B,, the large positive value, 459 A3,determined for Ar from the solubility measurements at 0.2 "C is consistent with the value, -13 A3, predicted for 25 OC. The theoretical calculations discussed in the body of the paper suggest that the nonideality of an aqueous solution of Kr in water is dominated by the second osmotic virial coefficient contribution even at concentrations up to about 1 M, which corresponds to mole fractions of order of 0.02. If this were true also for the similar
-
A3
(35) I. Kritchevsky and A. Iliinskaya, Acra Physicorhim. URSS, 20, 327 (1945). (36) W. L. Masterton, J . Chem. Phys., 22, 1830 (1954). (37) A. Ben-Naim, "Water and Aqueous Solutions", Plenum Press, New York, 1974. (38) "Handbook of Chemistry and Physics", 52nd ed., Chemical Rubber Co., Boca Raton, FL, 1971. (39) W. E. Deming and L. E. Shupe, Phys. Rev.,37, 638 (1931). (40) W. E. Deming and L. E. Shupe, Phys. Rev., 40, 848 (1932). (41) G . S . Kell, G. E. McLaurin, and E. Whalley, J . Chem. Phys., 68, 2199 I1 978) -... -,(42) D. R. Douslin, R. H. Harrison, R. T. Moore, and J. P. McCullough, J . Chem. Eng. Data, 9, 358 (1964). (43) A. Michels, Hub. Wijker, and Hk. Wijker, Physira, 15, 627 (1949). \ - -
802
J . Phys. Chem. 1986, 90, 802-807
solutes in Figure 6 , the plots should indeed be linear, as they are for all the solutes with the exception of N2, provided the data at the lowest mole fraction (approximately 1 X are discarded. Moreover, although the solubility data for Kr are not available, the correlation in Figure 7 strongly suggests that the larger and more soluble nonpolar solutes will have a positive osmotic virial
coefficient in water, as we have found in our theoretical calculations on Kr. This is experimental evidence for the assertion that some hydrophobic solutes do not experience hydrophobic attraction but rather a tendency to stay away from each other in solution. Registry No. Kr, 7439-90-9;water, 7732-1 8-5.
Excited-State Dynamics of Jet-Cooled Pyrene and Some Molecular Complexes Elisa A. Mangle and Michael R. Topp* Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104 (Received: August 8, 198s)
The fluorescenceexcitation spectra and fluorescencelifetimes for different levels of excitation have been measured for jet-cooled pyrene. The fluorescence lifetime decreases steadily with increasing energy from 1400 ns at the SI origin to -200 ns above 6500 cm-’ total vibrational energy. It is shown that the fluorescence lifetimes of S2 levels follow on continuously from the S1 levels, consistently with the strong vibronic coupling. However, the persistence of significant scatter in the fluorescence lifetimes indicates that the uncomplexed pyrene molecule is not yet in the statistical limit at an excess energy of 4000 cm-’. The fluorescence lifetime for pyrene complexes excited into the S I state is found to be the same or longer than that for corresponding excitation of the bare molecule. This shows that both directly excited high symmetry modes and the optically inactive modes reached by internal randomization of vibrational energy (IVR) of pyrene complexes are coupled to about the same degree to lower electronic states. This is quite different from the observed behavior of perylene complexes, where IVR causes a substantial reduction in fluorescence lifetime. The excitation spectrum of 1:1 n-pentane-pyrene complexes has the expected red shift on the S2resonances. However, a blue shift observed on the SI transitions for 1:l and 2:l complexes, evolving into a red shift at higher levels of aggregation, shows the evolution from molecular complexation, dominated by bimolecular interaction potentials, to larger aggregates, where an evolution toward a more crystalline type of environment appears to be taking place. Finally, it is shown that, because vibrational predissociation of the complexes generates long-lived fluorescence of free pyrene following excitation above the origin of Sz, clear fluorescence excitation spectra can be obtained for pyrene complexes in an otherwise congested region of the spectrum.
1. Introduction Recent work in this Iaborat~ryl-~ has studied the spectroscopic properties of isolated molecular complexes of the perylene molecule, generated in a supersonic free-jet expansion. In all cases studied, including alkanes, polyenes, alkyl and aryl halides, and benzene, red shifts have been observed of the complex resonances in relation to the corresponding bare molecule transitions. The red shift values measured have varied from
