Molecular Dynamics Simulation of Solvation Dynamics in Methanol

different solute sizes in mixtures with methanol mole fractions, xm = 0.2, 0.5, and 0.8, ... Inorganic Chemistry 2018 57 (16), 10050-10058 ... Loc...
4 downloads 0 Views 437KB Size
Download PDF
18258

J. Phys. Chem. 1996, 100, 18258-18268

Molecular Dynamics Simulation of Solvation Dynamics in Methanol-Water Mixtures Munir S. Skaf*,† Departamento de Quı´mica-FFCLRP, UniVersidade de Sa˜ o Paulo, AV. Bandeirantes, 3900 Ribeira˜ o Preto, SP 14040-901, Brazil

Branka M. Ladanyi* Department of Chemistry, Colorado State UniVersity, Fort Collins, Colorado 80523 ReceiVed: June 4, 1996; In Final Form: August 23, 1996X

The solvation dynamics following charge-transfer electronic excitation of diatomic solutes immersed in methanol-water mixtures is investigated through molecular dynamics simulations. The solvation response functions associated with an instantaneous reversal of the solute’s dipole moment for two different solute sizes in mixtures with methanol mole fractions, xm ) 0.2, 0.5, and 0.8, are calculated and compared to the corresponding ones in the pure liquids. The solvation response of the mixtures is separated into methanol and water contributions in order to elucidate the role played by each molecular species on the solvation dynamics. We find significantly different responses for the two solutes and relate them to the fact that the solute with the smaller site diameters is a much better hydrogen (H)-bond acceptor than the larger diameter solute. For the small solute in methanol and in the mixed solvents, we have also calculated H-bond response functions, which measure the rate of solute-solvent H-bond formation after the solute’s excitation and find that, at longer times, the solvation and H-bond formation response functions decay at similar rates. The implications of this finding for solvation dynamics of H-bonding solutes in H-bonding solvents are discussed and related to recent experimental results for such systems.

I. Introduction The dynamics of charge-transfer chemical reactions are often strongly influenced by the rate of solvent reorganization in response to the change in the solute charge distribution.1-3 There is therefore a considerable amount of interest in understanding the molecular mechanisms of solvation dynamics, i.e., of solvent relaxation in response to a sudden change in the solute charge distribution. This process can be detected in the absence of a chemical reaction by measuring the time-dependent Stokes shift of the fluorescence spectrum of a dissolved chromophore. There has been considerable recent progress4-7 in the measurement and modeling of solvation dynamics, stimulated in part by the advances in ultrafast spectroscopy, which have made it possible to detect this process with a time resolution of a few tens of femtoseconds. Molecular dynamics (MD) simulation has played an important role in predicting solvation response and determining its molecular origin. MD simulations show that the typical solvation response function, S(t), which measures the time evolution of the frequency at the maximum of the fluorescence band, exhibits an overall nonexponential decay in polar media.5,6 Its time evolution can be roughly cast into three regimes or stages of relaxation, namely, (1) an initial Gaussian (i.e., inertial) relaxation regime; (2) a subsequent librational regime, characterized by fast, hindered, oscillatory molecular motions, and (3) a slower diffusive regime. The relative importance of these three stages of solvation dynamics depends on solvent properties and on solute-solvent interactions before and after electronic excitation. A significant difference between hydrogen (H)bonding and nonprotic polar solvents is found in the onset and relative importance of the librational regime. Computer simula† Present address: Departamento de Fisico-Quimica, Instituto de Quimica, UNICAMP, Cx. Postal 6154, 13083-907 Campinas, SP, Brazil. X Abstract published in AdVance ACS Abstracts, November 1, 1996.

S0022-3654(96)01634-6 CCC: $12.00

tions of solvation dynamics in water9-16 and methanol,16-22 as well as our preliminary results for methanol-water mixtures,23 show that the duration of the inertial portion of S(t) is very short, around 20-40 fs, about 1/10 of its time scale in simple nonprotic polar solvents like acetonitrile22,24,25 and methyl chloride.26 The librational regime is much more prominent in H-bonding than in nonprotic solvents, with several periods of high-frequency oscillation, characteristic of the O-H libration subject to the H-bond restoring forces,27 appearing in S(t). Libration and other oscillatory hindered solvent motions also play an important role in solvation dynamics in complex liquids such as polyethers.28 Despite common characteristics, the solvation responses in pure water and methanol differ markedly from each other. The inertial and librational features of the water’s response are much more prominent and are responsible for most of the decay of S(t), while in methanol the diffusive regime plays an important role. The difference in the relative importance of inertial decay can be ascribed to the fact that all of water’s moments of inertia are small, while in methanol only one out of three moments of inertia is small.19,27,29 The differences in the librational and rotational-diffusion regimes of these two liquids are likely to be related to the difference in the number of H-bonds, as well as in the connectivity, strength, and persistence of their H-bond networks. Upon mixing, water and methanol exhibit interesting changes in structural, dynamical, and H-bond network properties.29-33 These changes and their molecular origin have been the subject of several computer simulation studies.27,29-33 Especially relevant to solvation dynamics is the composition-dependence of their dielectric properties.29,31,34 Our MD study of wave vector-dependent dielectric relaxation in methanol-water mixtures29 has revealed that these dynamics do not have a simple composition dependence due, in large part, to strong interspecies cross-correlations. In a recent article23 we have reported preliminary MD simulation results of solvation dynamics in an © 1996 American Chemical Society

Solvation Dynamics in Methanol-Water Mixtures

J. Phys. Chem., Vol. 100, No. 46, 1996 18259

equimolar methanol-water mixture using different sizes of diatomic probe solute molecules undergoing an instantaneous reversal (flip) of their dipole moments, (e/2,-e/2) f (-e/2,e/ 2). By separating the total response functions of the mixture into methanol and water contributions, we were able to identify methanol- and water-like features of the response. We found large differences in the equilibrium solvation structures for the two probe molecules, with the small probe exhibiting a much stronger H-bond acceptor tendency, and were able to relate them to the effects of the solute size on S(t). In the present work, we report a more thorough investigation of solvation dynamics in methanol-water mixtures, extending our studies to composition dependence by including methanol- and water-rich mixtures with “large” and “small” solutes. We also compare the mixture solvation dynamics results with the corresponding ones in the pure liquids. Further, we explore the effects of solute-solvent H-bonding on solvation dynamics by investigating the composition dependence of the solvation structure and by calculating the solute-solvent H-bond response functions in addition to S(t). It should be noted that several experimental studies of solvation dynamics in alcohols point to the contributions of solute-solvent H-bonding to this process.35-38 For example, Jonah and co-workers35 found that solvation dynamics of the benzophenone anion in various alcohols is considerably faster than predicted on the basis of nonspecific solvation models.5,7,8 Evidence of specific solvation due to solute-solvent H-bond formation has been found by Berg and co-workers36 for the absorption recovery of resorufin in ethanol and ethylene glycol and by Rullie`re and co-workers35 for the fluorescence of the chromophore MPQB in 1-propanol. Maroncelli and co-workers37 studied the solute dependence of solvation dynamics in 1-propanol and found that S(t) of some classes of chromophores differs considerably from what one would expect on the basis of nonspecific solvation models and that this difference can be correlated with the likely magnitude of the change in solutesolvent H-bond energies on electronic excitation. We hope that our study, which allows us to analyze the effects of H-bonding in greater detail than is possible in experimental investigations, will help elucidate this important aspect of solvation in protic polar solvents. The remainder of this article is organized as follows. In the next section, we present the essential details of the simulations and potential parameters, define the solvation response function separated into methanol and water contributions, and introduce the solute-solvent H-bond formation response function. In section III, we present our results. Our main conclusions are summarized in section IV. II. Response Functions and Computational Details A. Solvation and H-Bond Response Functions. In ultrafast time-dependent fluorescence Stokes shift (TDFSS) experiments, one measures the temporal evolution of the emission frequency from a charge-transfer excited electronic state of a dissolved chromophore. The Stokes shift response is usually reported as the normalized response function of the frequency, ν(t), at the maximum of the fluorescence band4-7

S(t) )

hν(t) - hν(∞) hν(0) - hν(∞)

(1)

In this expression t ) 0 corresponds to the time of electronic excitation and t ) ∞ corresponds to the time at which steadystate Stokes shift has been reached. Since the electronic excitation occurs essentially instantaneously on the time scale of nuclear motions, the time evolution of the Stokes shift represents the dynamics associated with rearrangements of the

nuclei in response to the charge redistribution within the probe. If the probe molecule is, to a good approximation, rigid, these nuclear rearrangements are due to solvation dynamics, i.e., to the time evolution of the difference, ∆E, between solute-solvent interaction energies for the ground (G) and excited (E) electronic states.39 In our simulations the solute is a diatomic, and its electronic excitation is represented by an internal redistribution of the solute atomic charges leading to a sudden reversal of the solute’s dipole moment: (e/2,-e/2) f (-e/2, e/2). The energy gap is given then by the electrostatic solute-solvent interactions:

∆E ) UelE - UelG

(2)

where UEel and UGel are the excited- and ground-state solutesolvent electrostatic potential energies, which are readily obtained from the simulations. The solvation response function is given by38

S(t) )

∆E(t) - ∆E(∞) ∆E(0) - ∆E(∞)

(3)

where the overbar indicates an average over nonequilibrium trajectories, which are generated by propagating the system in the presence of the excited-state solute, starting from several independent equilibrium configurations with the solute in its ground state. Due to the symmetry between the ground and the excited solute states, the average ∆E(∞) can be calculated from the equilibrium trajectories without having to specify the state the solute is in. In order to identify methanol- and waterlike features to the total solvation response, we have calculated separately the solute-methanol, ∆Em, and solute-water, ∆Ew, contributions to the interaction energy difference ∆E. The normalized response function is then separated into methanol and water contributions:

S(t) ) Sm(t) + Sw(t)

(4)

where

SR(t) )

∆ER(t) - ∆ER(∞) ∆E(0) - ∆E(∞)

(5)

and R ) m, w. In addition to the solvation response functions, we also study the rate of solute-solvent H-bond formation after the solute’s excitation. For this purpose we introduce the following normalized H-bond response function:

SHB(t) )

nHB(t) - nHB(∞) nHB(0) - nHB(∞)

(6)

where nHB(t) is the average number of solute-solvent H-bonds at time t calculated from the nonequilibrium trajectories. Again, due to the symmetry between the ground and the excited solute states, there is no need to specify the solute’s state in order to calculate the average nHB(∞) from the equilibrium configurations. As in the case of S(t), one can calculate separate contributions from solute-methanol and solute-water H-bonds to the total response SHB(t), namely, m w (t) + SHB (t) SHB(t) ) SHB

where

(7)

18260 J. Phys. Chem., Vol. 100, No. 46, 1996

S

R

HB(t)

)

Skaf and Ladanyi

R R nHB (t) - nHB (∞)

(8)

nHB(0) - nHB(∞)

TABLE 1: Site Parameters for the Diatomic Solutes Used in the Simulations solute

R/kB/K

σR, Å

qR, e

mR, au

small large

87.9 87.9

3.083 4.200

(0.5 (0.5

30 30

with

a

m w (t) + nHB (t) nHB(t) ) nHB

(9)

m w (t) and nHB (t) are respectively the number of The quantities nHB solute-methanol and solute-water H-bonds at time t. These can be calculated from the trajectory data using a geometric definition for H-bonding27,30 between the negative solute site and the solvent hydroxyl hydrogen, as described in part B of this section. B. Interaction Potentials and Simulation Details. In our simulations we have studied probe diatomic solute molecules19,23 with atomic partial charges e/2 and -e/2 immersed in methanolwater mixtures of different compositions as well as in pure water. The molecules interact via site-site Lennard-Jones plus Coulomb potentials, obeying Lorentz-Berthelot combining rules, a common way of approximating interactions between unlike sites on different molecules.40 Specifically, the potential between sites R and β on different molecules is

VRβ(r) ) 4(Rβ)1/2

[(

) (

σR + σβ 2r

1/2

-

)]

σR + σβ 2r

6

+

qRqβ 4π0r (10)

where r is the site-site distance, qR is the partial charge on site R, and R and σR are the Lennard-Jones (LJ) potential well depth and diameter. Both solute sites have the same values of R and σR. We consider solvation dynamics of two solute diatoms which differ from each other only by the values of the LJ diameter σR. We refer to the two solutes as “small” (σR ) 3.083 Å) and “large” (σR ) 4.20 Å). The small solute site diameters are the same as those of methanol oxygens, while the large site diameters exceed the size of the methyl group and would also exceed the size of large atomic sites such as iodine, typically set at about 4 Å.41 As we shall see, the most significant consequence of this size difference is in its effect on the H-bonding ability of the two solutes. The values of the solute site parameters, R, σR, qR, and masses, mR, are given in Table 1. The solute molecular bond length is fixed at 1.40 Å, slightly different from the 1.43 Å value used in refs 18 and 19. For methanol we use the H1 model42 and for water the TIP4P model.43 The site parameters and molecular geometries for both models may be found in ref 30. In both models, the hydrogens of the O-H group carry partial charges but do not have LJ potentials associated with them. The TIP4P model has a single LJ potential site coincident with oxygen, while the H1 model has two LJ sites: one corresponding to O and the other to the methyl group, which is represented as a single site. The TIP4P and the H1 models have both been thoroughly tested and are known to provide good representations of the thermodynamic, structural, and transport properties of the pure liquids near room temperature.27,42,43 The models also predict static dielectric properties and far-infrared absorption spectra in good agreement with experiment,44-46 while underestimating somewhat the Debye relaxation times.45,46 Mixtures of TIP4P water and H1 methanol near room temperature have been extensively studied as well.29-31 Good agreement with experiment has been obtained for the excess thermodynamic quantities of mixing, diffusion coefficients, and H-bond lifetimes.30 The agreement with experiment for the dielectric properties is comparable in quality to that for the pure components.29,31 For methanolwater mixtures experimental structural data are unavailable, but

The bond length is 1.40 Å for both solute types.

TABLE 2: Simulation Detailsa xm

Nm

Nw

box length L, Å

density, Å-3

no. of trajectories

0.0 0.2 0.5 0.8

0 100 250 399

255 399 249 100

19.73 26.30 28.66 30.92

0.03450 0.02749 0.02124 0.01691

200 200 240 240

a

These values correspond to both solute types.

MD results, based on potentials of the type given by eq 10, for other aqueous mixtures with molecules capable of acting as acceptors (DMSO47) or as both donors and acceptors (ammonia30) of H-bonds have been successfully compared with neutron and X-ray scattering data. All this gives us confidence that the present model for methanol-water mixtures is reasonable and that it will provide us with insight into solvation dynamics in such mixtures. For each solute size, three sets of MD simulations in the microcanonical ensemble with periodic boundary conditions were carried out on systems of 500 molecules (one of them being the diatomic solute probe) corresponding to a methanolrich mixture, with methanol mole fraction xm ) 0.8, an equimolar mixture, and a water-rich mixture, with xm ) 0.2, at an average temperature of 298 K. In addition, we also simulated solvation dynamics in pure TIP4P water, using systems of 1 solute and 255 water molecules. The sizes of the cubic simulation boxes were determined from the experimental densities of the mixtures at atmospheric pressure and are given in Table 2. The cutoffs for the Lennard-Jones interactions were set to half the box length for each mixture, and the long-range part of the Coulombic forces was calculated using Ewald sums39 with conducting boundary conditions. The equations of motion were integrated using the leapfrog algorithm39,48 with a time step of 4 fs for the mixtures and 3 fs for pure water. The molecular geometries of the solute molecule and of methanol were maintained using SHAKE,49 and that of water was maintained by the generalized method of constraints.50 As shown in ref 50, constraint methods are more stable than quaternion-based MD algorithms for rigid molecules, allowing larger time steps for a given level of accuracy. In our equilibrium MD simulations, we find that for the equimolar mixture the relative energy fluctuations are of (0.017% in an 8 ps run, and for pure water, we find that the energy fluctuations are (0.014% in a 6 ps run. The equilibrium data that we present for a given solute and mixture composition are averages of 9-11 trajectories of 8000 time steps, separated by 800 time step constant-temperature runs. The nonequilibrium simulations were performed as follows: 18-20 first, a set of statistically independent initial configurations was extracted from the equilibrium runs with the solute in its ground state; then, an instantaneous switch of the solute’s atomic charges was performed, and the system was allowed to relax to equilibrium. During this process, the energy gaps ∆Em and ∆Ew, as well as the molecular coordinates were recorded for later analysis. At least 200 independent nonequilibrium trajectories of 1.2 ps duration were recorded for each solute type and mixture and used to calculate the average ∆E and nHB values entering eqs 4, 5, 7, and 8. Additional details about the simulations are given in Table 2. In order to identify solute-solvent H-bonds, we use geometrical criteria analogous to those used for water,27 methanol,42

Solvation Dynamics in Methanol-Water Mixtures

J. Phys. Chem., Vol. 100, No. 46, 1996 18261

Figure 1. Solute-water (dashed line) and solute-methanol (full line) pair distributions g-HR(r) and g-OR(r) involving the negative site of the small solute in the equimolar methanol-water mixture. g-HR(r) is shifted vertically by 3. The geometric criteria identifying solute-solvent H-bonds are shown schematically.

and their mixtures.30 They are based on the locations of the first peaks in pair distribution functions g-HR(r) and g-OR(r), where - denotes the negative solute site. These functions for the small solute in the equimolar methanol-water mixture are shown in Figure 1. Since the solute sites are similar in size to the solvent oxygen sites, we use the same distance criteria as for O-H- - -O hydrogen bonds,30 specifically r-HR < 2.6 Å, r-OR < 3.5 Å and the angle ∠-ORHR < 30°. As can be seen from Figure 1, the criterion corresponds to the locations of sharp first peaks in g-HR(r) and g-OR(r) distribution functions. As we shall demonstrate in our more detailed presentation of the results on the solvation structure, g-HR(r) for the large solute is quite different. It exhibits a much smaller and more diffuse first peak, indicating that this solute is a much weaker H-bond acceptor. III. Results and Discussion We start by describing our results for S(t) and its components, Sm(t) and Sw(t), for methanol-water mixtures of different composition. We display next the results for the solvation structures around the two solute types and then compare the solvation and solute-solvent H-bond formation dynamics. A. Solvation Responses. In Figure 2 we display the results for S(t) for dipole flip in the small and large solutes. In addition to our results for xm ) 0.8, 0.5, 0.2, and 0.0, we show the pure methanol data, which were previously obtained by Fonseca and Ladanyi.19 For both solute types, S(t) is highly nonexponential. Its initial Gaussian decay portion increases in importance and decreases in duration as the concentration of water increases. This trend can be understood by noting that all of water’s moments of inertia are small, corresponding predominantly to hydrogen motion. In methanol, only one out of three moments of inertia is similarly small.27,29 A water molecule can also participate in a larger number of H-bonds than methanol.30 This affects the librational portion of S(t). As the concentration of water increases, the high-frequency oscillatory features, characteristic of H-bond librations, become more pronounced. The librational and diffusive portions of S(t) are markedly different for the two solutes. In the case of the large solute, the amplitude of librational oscillations is significantly more pronounced and the subsequent diffusive decay more rapid than for the small solute. The small solute response becomes gradually slower

Figure 2. Solvation response functions S(t) for dipole flip, (e/2,-e/2) f (-e/2,e/2), of the large (top panel) and small (bottom panel) solutes in pure methanol (thin solid line, from ref 19), pure water (thick solid line), and methanol-water mixtures: xm ) 0.8 (dashed line), xm ) 0.5 (dash-dotted line), and xm ) 0.2 (dash-double dotted line).

as xm increases, while S(t) for the large solute exhibits a more complicated composition dependence. For t > 0.1 ps, S(t) in pure methanol decays more rapidly than in mixtures with xm ) 0.8 and 0.5. For t > 0.25 ps S(t) in the equimolar mixture becomes more slowly decaying than in the methanol-rich mixture. These trends parallel the ones we observed for the composition dependence of the longitudinal dipole densitytime correlations in methanol-water mixtures.29 This suggests that the large solute does not strongly perturb the solvent structure, and its solvation dynamics reflects collective dipolar fluctuations in the pure solvent. These, in turn, are strongly influenced by the changes with xm in the concentration and connectivity of the H-bond network in methanol-water mixtures. A clearer picture of the solute dependence of solvation dynamics is shown in Figure 3, where the solvation responses for the two solutes are compared directly. These data show that the difference in the solvent responses to the two solutes diminishes as the concentration of water increases. The most significant and pronounced differences occur in the librational portion of the response. In the case of the large solute, the onset of H-bond libration occurs sooner and lasts longer than for the small solute. These differences can be explained by noting that the initial solvation structures are quite different for the two solutes. As we have shown for the equimolar mixture, the small solute has a much stronger tendency to form a hydrogen bond with the solvent. The initial step in its solvation dynamics is the rapid destruction of this H-bond when the solute site charges are reversed. As Fonseca and Ladanyi18,19 have shown for solvation in pure methanol, this occurs mainly through the free rotation of the affected hydrogen around the appropriate axis (e.g., the C-O bond in methanol) of the solvent molecule. Since the large solute has a much smaller tendency to form

18262 J. Phys. Chem., Vol. 100, No. 46, 1996

Skaf and Ladanyi

Figure 3. Comparison of the solvation responses for dipole flip of the large (solid line) and small (dashed line) solutes. The methanol mole fractions are as indicated. The responses for xm > 0 are shifted vertically by multiples of 0.3.

H-bonds with the solvent, this free rotation step plays a greatly diminished role in its solvation dynamics. The data in Figure 3 also indicate that the diffusive portion of S(t) is slower for the small solute. As we shall see, much of this portion of the response corresponds in this case to the formation of a new solute-solvent H-bond. The data in Figure 3 show that the solvation responses of the two solutes become increasingly similar as the concentration of water is increased. This increase in similarity of the responses also indicates that a better agreement with linear response theory is achieved as xw increases. Fonseca and Ladanyi had shown that linear response provides a poor approximation for dipole flip19 and dipole creation18,19 in pure methanol, especially for the small solute for which these changes represent larger perturbations in the solvation energy. The best agreement with linear response can be expected in pure water, based on MD results for ionic solvation dynamics in this solvent.9,11 Figure 4 illustrates the comparison of S(t) for dipole flip in pure TIP4P water with its linear response counterpart

C(t) ) 〈δ∆E(0) δ∆E(t)〉/〈(δ∆E)2〉

(11)

where δ∆E ) ∆E - 〈∆E〉 denotes a fluctuation in the solvation energy difference between excited and ground state solutes and 〈...〉 indicates an equilibrium ensemble average. Note that since the solute-solvent potentials are the same for the ground and excited state solutes, the solute state for which C(t) is evaluated need not be specified. As can be seen from Figure 4, C(t) and S(t) are far from being in perfect agreement, even in pure water. For both solutes, the initial decay of C(t) is slower than that of S(t) and the inertial contribution to C(t) is considerably smaller, around 40% vs over 60% for S(t). It thus appears free rotation of water molecules, which is responsible for the rapid initial decay of S(t), plays a less important role in the short-time dynamics of C(t). This seems reasonable in view of the expectation that most hydrogens in the vicinity of the solute are likely to participate in H-bonds and at equilibrium would not be free to rotate. Free rotation

Figure 4. Comparison between the nonequilibrium solvation response S(t) (solid line) and the δ∆E time correlation, C(t) (dashed line), its linear response counterpart, for dipole flip in pure water. The results for the large and small solutes are shown respectively in the upper and lower panels.

can occur after the solute dipole reversal leads to breakup of some of H-bonds. As noted earlier, this nonequilibrium solvation mechanism is more important for the small solute, given that it is initially more likely to be H-bonded to the solvent. We turn now to the analysis of S(t) in terms of the relative contributions from the two cosolvents. Before discussing the analysis of S(t), we discuss the effects of mixing on the average solvation energy differences. This information is given in Table 3 which contains values of 〈∆ER〉, R ) m, w, for the three mixtures and for pure water and methanol. In the present case of the dipole flip,

∆ER(0) - ∆ER(∞) ) -2〈∆ER〉

(12)

Table 3 shows that 〈∆Ew〉 in pure water is more attractive, i.e., a larger negative number, for both solute types than 〈∆Em〉 in pure methanol. Denoting these values by subscript p (for pure), we find

〈∆Ew〉p ) 1.29〈∆Em〉p

(13)

As proposed by Chandra and Bagchi in their work on solvation in dipolar mixtures,51 the extent of nonideality in equilibrium solvation can be estimated by calculating how much 〈∆Ew〉 and 〈∆Em〉 in mixtures deviate from what would be expected on the basis of composition and pure component average solutesolvent electrostatic energies. Specifically, this deviation, which we call the excess solvation energy difference, is given by

∆〈∆ER〉ex ) 〈∆ER〉 - xR〈∆ER〉p In Table 3, we list the values of ∆〈∆ER〉ex as well as of

(14)

Solvation Dynamics in Methanol-Water Mixtures

J. Phys. Chem., Vol. 100, No. 46, 1996 18263

TABLE 3: Average Equilibrium Values of the Electrostatic Solute-Solvent Energy Gap small solute xm 0.0 0.2 0.5 0.8 1.0

〈∆Em〉, kJ/mol -21.7 -53.7 -77.6 -100.4

∆〈∆Em〉ex, kJ/mol -1.6 -3.5 2.7 0.0

∆ m, % 8.0 7.0 -3.4 0.0

∆R ) 100

large solute

〈∆Ew〉, kJ/mol

∆〈∆Ew〉ex, kJ/mol

∆w, %

-129.6 -104.0 -59.4 -24.7

0.0 -0.3 5.4 1.2

0.0 0.3 -8.3 -4.6

∆〈∆ER〉ex xR〈∆ER〉p

〈∆Em〉, kJ/mol -13.4 -29.4 -33.3 -42.5

∆〈∆Em〉ex, kJ/mol -4.9 -8.2 0.7 0.0

∆m, % 58 39 -2.1 0.0

〈∆Ew〉, kJ/mol

∆〈∆Ew〉ex, kJ/mol

∆w, %

-54.9 -39.7 -18.8 -11.6

0.0 4.2 8.7 -0.6

0.0 -9.6 -32 5.5

(15)

which measures the relative excess solvation energy difference. These data indicate that both quantities vary with composition, solute, and solvent type. For both solutes, ∆〈∆ER〉ex values are the largest for the equimolar mixture, where the interaction with methanol is more attractive than expected on the basis of eq 13 and the interaction with water less attractive. Larger deviation in both absolute and relative terms is found for the large solute. This solute also shows a dramatically enhanced affinity for methanol at xm ) 0.2, accompanied by a suppressed affinity for water. At xm ) 0.8 the large solute exhibits a slight preference for water over methanol. In the case of the small solute at xm ) 0.2 and 0.8, the interactions with the two solvents do not compensate for each other. Thus, at xm ) 0.2 both ∆〈∆ER〉ex are negative, indicating that interactions with both solvents are more attractive than expected, while ∆〈∆ER〉ex is positive at xm ) 0.8, showing that they are less attractive for both water and methanol. Overall, the small solute behaves more ideally than the large solute, which shows a strong preference for methanol in equimolar and water-rich mixtures. In Figures 5 and 6 we display S(t), Sm(t), and Sw(t) for the three mixtures for the large and small solutes, respectively. The initial, t ) 0, values of Sm(t) and Sw(t) are already anticipated from the data shown in Table 3. At longer times the most dramatic difference in solvation dynamics of the two solutes occurs in the equimolar mixture, where the contributions of methanol and water are approximately equal in the case of the small solute, while Sm(t) remains appreciably larger than Sw(t) for the large solute over the entire time interval shown. For the large solute, methanol’s contribution to S(t) is almost as important in the equimolar as in the methanol-rich mixture. By contrast, despite the preference of this solute for methanol in the water-rich mixture, the longer time decay of S(t) in this mixture is strongly dominated by water. This appears to be due to the rapid decay of Sm(t) in this mixture. In Figures 7 and 8 we examine the composition dependence of the decay rates of Sm(t) and Sw(t) for the two solutes. Figure 7 shows that Sw(t) for the large solute changes gradually in going from pure water to the equimolar mixture, with the slight decrease in the initial curvature and in the inertial component amplitude. In going from xm ) 0.5 to xm ) 0.8 both the curvature and the inertial component amplitude decrease much more dramatically, indicating that the water’s motion has become much more strongly hindered in the methanol-rich mixture. Figure 7 shows that Sm(t) of pure methanol and in the xm ) 0.8 are essentially identical. In the equimolar mixture, the amplitude of the inertial component increases slightly, while the diffusive decay becomes slower than at xm > 0.5. Note that the long-time decay of S(t) for the large solute exhibits the same trends as Sm(t) for xm ) 0.5 and 0.8. Sm(t) in the waterrich mixture is very different from the methanol responses at other compositions. It lacks any librational features, and it decays considerably more rapidly at longer times. This and the fact that it decays more rapidly than Sw(t) at xm ) 0.2 suggests that methanol molecules in the vicinity of the solute

Figure 5. Total solvation response functions (solid lines) and their contributions from methanol, Sm(t) (dashed lines), and water, Sw(t) (dashdotted lines), for the large solute in methanol-water mixtures.

in this mixture may not be efficiently incorporated into the solvent H-bond network. Figure 8 shows that Sw(t) and Sm(t) exhibit less dramatic composition dependence than in the case of the large solute. Both the water and methanol responses show small decreases in the initial curvatures, inertial component amplitudes, and diffusive decay rates as the methanol concentration increases. Comparing the methanol and water responses shown in Figures 7 and 8, we see that librational oscillations are considerably more pronounced for Sw(t) than in Sm(t), as we would expect in view of the larger number of H-bonds that water makes. Addition of methanol leads to only a slight damping of librational oscillations in Sw(t). As noted earlier for pure methanol,19 Sm(t) for the small solute lacks any oscillatory librational features. For the large solute these features are absent only in the water-rich mixture. B. Solvation Structures. We have seen that there is a considerable difference in the way that S(t) varies with composition for the two solutes and that both Sw(t) and Sm(t) exhibit more dramatic composition dependence for the large than for the small solute. Additional insight into the solvation of the two solute types can be gained by examining the solvent structure around them. Of the 14 solute-solvent site-site pair

18264 J. Phys. Chem., Vol. 100, No. 46, 1996

Figure 6. Same as Figure 5, but for the small solute in methanolwater mixtures.

distributions for each solute type, we have chosen to display only six: g-Hw(r), g-Hm(r), g+Ow(r), g+Om(r), g-Me(r), and g+Me(r). The first four correspond to structures that are governed primarily by attractive solute-solvent electrostatic interactions, while the remaining two correspond to structures that are primarily nonelectrostatic in nature, given that the methyl group is the largest solvent site carrying a large LennardJones energy (Me/kB ) 91.2 K) and a weak positive charge (qMe ) 0.297e). Because both solvents contain O and H sites, we have grouped on the same figures pair distributions of a solute site with these similar sites on the two solvents. Thus, Figures 9 and 10 both display g-Hw(r) and g-Hm(r), the first one for the large solute and the second for the small one. Similarly, Figures 11 and 12 include g+Ow(r) and g+Om(r) for the large and small solutes, respectively. Finally, Figures 13 and 14 depict respectively g-Me(r) and g+Me(r) for both solutes. In comparing Figures 9 and 10, we see that for a given solute g-Hw(r) and g-Hm(r) are quite similar but that they are strikingly different for the two solutes. The difference stems from the fact that the small solute is a much better H-bond acceptor than the large solute. Thus, g-Hw(r) and g-Hm(r) for the small solute exhibit a tall and sharp first peak centered at around 2 Å, while the first peak in the large solute pair distributions at around 2.3 Å is much smaller and more diffuse. Figure 9 shows that, at all compositions, the second peak in g-Hw(r) and g-Hm(r) is taller than the first, indicating that the hydroxyl hydrogen tends to be oriented away from the negative site of the large solute. g-Hm(r) for the large solute varies gradually with xm, with the second peak moving inward and increasing in height as xm decreases. By contrast, both the second and third peaks in g-Hw(r) have a complicated composition dependence. For example, the second peak in the water-rich mixture is quite different from the pure water result.

Skaf and Ladanyi

Figure 7. Normalized water (upper panel) and methanol (lower panel) solvation responses for the large solute. Sw(t)/Sw(0) is shown for xm ) 0.0, 0.2, 0.5, and 0.8; Sm(t)/Sm(0) is shown for xm ) 0.2, 0.5, 0.8, and 1.0.

Figure 8. Same as Figure 7, but for the small solute.

In the case of the small solute, the peak locations in g-Hw(r) and g-Hm(r) are weakly composition dependent. The heights

Solvation Dynamics in Methanol-Water Mixtures

J. Phys. Chem., Vol. 100, No. 46, 1996 18265

Figure 9. Pair distributions involving the negative site (-) of the large solute and the hydrogen (H) solvent sites. g-Hw(r) for xm ) 0.0, 0.2, 0.5, and 0.8 is depicted in the bottom panel, and g-Hm(r) for xm ) 0.2, 0.5, 0.8, and 1.0 is depicted in the top panel.

Figure 11. Pair distributions involving the positive site (+) of the large solute and the oxygen (O) solvent sites. g+Ow(r) for xm ) 0.0, 0.2, 0.5, and 0.8 is depicted in the bottom panel and g+Om(r) for xm ) 0.2, 0.5, 0.8, and 1.0 is depicted in the top panel.

Figure 10. Same as Figure 9, but for the small solute.

Figure 12. Same as Figure 11, but for the small solute.

of the first peaks decrease as xm decreases. However, as we shall see, the number of solute-solvent H-bonds actually

increases with increasing water concentration. One should note that as xm goes from 0.8 to 0.2, the number of O-H bonds

18266 J. Phys. Chem., Vol. 100, No. 46, 1996

Skaf and Ladanyi TABLE 4: Average Number of Solute-Solvent Hydrogen Bonds for the Small Solute

Figure 13. Pair distributions involving the negative site (-) of the large (top panel) and small (bottom panel) solutes and the methyl group (Me) of methanol. g-Me(r) is shown for xm ) 0.2, 0.5, 0.8, and 1.0.

Figure 14. Same as Figure 13, but for the pair distribution g+Me(r) involving the positive solute site and the methyl group.

increases by a factor of 1.5. In the case of g-Hw(r), the heights of the second and third peaks also decrease with increasing water concentration. The third peak corresponds to an Hw that forms an H-bond to the oxygen of the water molecule H-bonded to

xm

m 〈nHB 〉

w 〈nHB 〉

〈nHB〉

0.2 0.5 0.8

0.11 0.41 0.62

1.25 0.67 0.27

1.36 1.08 0.89

the solute. The increase in its height with increasing xm can be attributed to the increased structuring with dilution of water in methanol-water mixtures.30 Using the geometrical criteria illustrated in Figure 1, we have calculated the number of solvent-solute H-bonds for the small solute in the three mixtures. The results are given in Table 4. A similar calculation for the large solute showed that 〈nHB〉, the average number of solute-solvent H-bonds, was very small, about 0.1 in all mixtures, so these results are not reported in detail. It can be seen from Table 4 that the number of solute-solvent H-bonds increases with increasing water concentration. Noting that water can supply two hydrogens and methanol only one, m w 〉/〈nHB 〉 ) 0.125, 0.50, and 2.0 in we would expect 〈nHB mixtures with xm ) 0.2, 0.5, and 0.8, based solely on the Hm and Hw concentrations. However, these ratios are actually 0.088, 0.609, and 2.3, showing that solute-methanol H-bonds are slightly favored in equimolar and methanol-rich mixtures, while solute-water H-bonds are slightly favored in water-rich mixtures. We should point out that we also find that, while the fluctuations in nHB are relatively small, significantly larger m w and nHB , indicating that the solute fluctuations occur in nHB does not have a strong preference for forming H-bonds with a particular solvent in methanol-water mixtures. We turn now to g+Ow(r) and g+Om(r), shown in Figure 11 for the large solute and in Figure 12 for the small solute. Again, more pronounced changes with composition occur for the large solute. For the small solute, g+Om(r) is weakly compositiondependent. The major trend is the decrease in its first peak and the increase in its second peak as xm decreases. For the large solute, g+Om(r) is quite unstructured for xm g 0.5 and develops a broad peak, centered at around 5 Å, for xm ) 0.2. This indicates that close proximity of Om to the positive solute site is not favored. By contrast, Ow is likely to be close to the positive solute sites of both large and small solutes, as the data on g+Ow(r) indicate. For both solutes, g+Ow(r) in mixtures has a sharp first peak whose height decreases with increasing water concentration. The peak is sharper and higher for the small solute, as one might expect based on its smaller site diameter and the resulting stronger attraction between its positive site and the negative charge site in water. We do not anticipate the solute interactions with the methyl group to be primarily electrostatic. Indeed, comparison of Figures 13 and 14 indicates that the structures of the methyl group around the negative solute site do not differ greatly from the structures around the positive site. The main difference is in the first peak, which is slightly narrower and higher in g-Me(r), which corresponds to a pair of oppositely charged sites, than in g+Me(r). The most striking aspect of Figures 13 and 14 is the difference in the composition dependence of the solutemethyl pair distributions for the two solute types. In the case of the small solute, g-Me(r) and g+Me(r) are nearly independent of composition, while for the large solute the height of their first peaks grows substantially as methanol concentration decreases from xm ) 1.0 to 0.2. This behavior of g-Me(r) and g+Me(r) indicates that the large solute is solvated predominately by the weakly polar methyl group. As water is added, the probability of finding a methyl group in the first solvation shell increases. However, at low methanol concentration, it becomes impossible to maintain a high concentration of methyl groups

Solvation Dynamics in Methanol-Water Mixtures

J. Phys. Chem., Vol. 100, No. 46, 1996 18267

Figure 16. Comparison of the decay rates of solvation and H-bond responses for the small solute in methanol-water mixtures. Shown are ln[S(t)] (solid line) and ln[SHB(t)] (dashed line). ln[SHB(t)] is shifted downward with respect to the corresponding ln[S(t)] by 1.1, 0.8, and 0.8 at to xm ) 0.8, 0.5, and 0.2.

Figure 15. Solute-solvent H-bond response functions for the small solute in methanol-water mixtures. Shown is the total H-bond response, m SHB(t) (solid line), and its contributions from methanol, SHB (t) (dashed w line), and water, SHB (t) (dash-dotted line).

in the first solvation shell. This explains why S(t) for the large solute changes drastically at xm e 0.2. In the case of the small solute, the composition of the first solvation shell is determined to a much larger extent by attractive electrostatic interactions, including formation of H-bonds with Hw and Hm. Thus, the composition of the first solvation shell does not differ greatly from the bulk composition of the solvent mixture, explaining the gradual change of S(t) with xm for this solute. C. Solute-Solvent H-Bond Responses. As we have shown, an important feature of the solvent structure around the small solute is the solute-solvent H-bond. In the case of dipole flip in methanol, Fonseca and Ladanyi19 showed that the destruction of this bond contributes to the fast initial decay of S(t) for this solute. They also showed, by examining the time evolution of g-Hm and g-Om, that the new solute-solvent H-bond forms slowly. Phelps et al.20 have identified in the same way H-bond formation as the slow step in anion-neutral charge migration in methanol. Here we investigate the H-bond formation process more directly by calculating the solute-solvent H-bond response functions, defined by eqs 7 and 8. The results for the three mixtures are displayed in Figure 15. The figure shows that the rate of solute-solvent H-bond formation increases as the concentration of water increases. In the water-rich mixture, SHB(t) is due almost entirely to w (t), as could be anticipated from Table 4 which indicates SHB that only about 8% of the solute-solvent H-bonds in this mixture at equilibrium are due to methanol. For the other two mixtures, both methanol and water make appreciable contribuR 〉 values tions, roughly in the proportions predicted by 〈nHB given in Table 4.

Comparison of Figures 6 and 15 indicates that at short times S(t) varies considerably faster than SHB(t). We do expect, however, that the H-bond formation process will make a significant contribution to the longer-time portion of S(t). We investigate this in Figure 16 by comparing the decay rates of S(t) and SHB(t), plotted on the logarithmic scale. We see that the longer-time decay rates of S(t) and SHB(t) are remarkably similar, suggesting that a substantial portion of the diffusive decay of S(t) can be ascribed to solute-solvent H-bond formation. IV. Summary and Conclusion In this article we have reported MD simulation results for solvation dynamics in methanol-water mixtures in response to an instantaneous reversal of the dipole moment of diatomic probe molecules. We have compared solvation responses of mixtures of different compositions (xm ) 0.2, 0.5, and 0.8) to solvation dynamics in pure water and pure methanol20 and analyzed it in terms of the relative contributions of the two cosolvents. We find very different solvation dynamics for the “small” and “large” solute diatoms, which differ only in the sizes of their site diameters. S(t) for the small solute changes gradually as water is added, relaxing faster as xm decreases. For this solute the component responses Sm(t) and Sw(t) also speed up gradually, but less than the total S(t), as water is added. By contrast, S(t) for the large solute shows a complicated composition dependence: Its longer-time decay rate is faster for pure methanol than for mixtures with xm ) 0.8 and 0.5, while its short-time (t < 40 fs) behavior changes drastically as xm changes from 0.5 to 0.2. Decomposition of S(t) for this solute indicates that it is dominated by Sm(t) for xm g 0.5. In order to account for the observed solute-dependence of S(t) and of its components, we have analyzed the composition dependence of the solvent structure around both solutes. We find that the small solute is a good H-bond acceptor which therefore has a similar affinity for both solvents. The solutesolvent pair distributions indicate that the large solute is

18268 J. Phys. Chem., Vol. 100, No. 46, 1996 hydrophobic and that it exhibits a similar lack of affinity for the hydroxyl hydrogen of methanol. To the extent allowed by the mixture composition, this solute is primarily solvated by the methyl group. This accounts for the drastic changes in S(t) at 0.2 when this type of solvation is no longer possible. For solvation in methanol, it has been shown earlier19,20 by looking at the time evolution of the solvent structure around the solute that solute-solvent H-bond formation is slow on the solvation time scale. We have here measured more directly the rate at which this process occurs for the small solute in methanolwater mixtures. By comparing the time evolution of logarithms of S(t) and SHB(t), we have shown that the diffusive decay rates for the two processes are essentially the same. Given also that S(t) decays more slowly in the post-librational regime for the H-bonding small solute than for the nonspecifically solvated large one, we conclude that solute-solvent H-bond formation is an important slow step in solvation dynamics of H-bonding solutes in protic polar solvents. In view of experimental results which have also identified solute-solvent H-bond formation as a source of solute dependence of solvation dynamics in H-bonding solvents,36,38 this process merits further investigation. MD simulation of solvation dynamics of model solutes which resemble more closely the chromophores used in experimental studies would be desirable. For example, a study focusing on the sources of differences in S(t) for aromatic amines and coumarins in alcohols would show to what extent the solutesolvent H-bond formation process that we have identified here for simple model solutes is relevant for experimentally accessible solvation probes. We are presently working on such simulations. Acknowledgment. This work has been supported by Grants CHE-9019707 and CHE-9520619 from the National Science Foundation and by funds from the Brazilian agency FAPESP through Grant 94/0273-4. References and Notes (1) Heitele, H. Angew. Chem., Int. Ed. Engl. 1993, 32, 359. (2) Hynes, J. T. In Ultrafast Dynamics of Chemical Systems, Simon, J. D., Ed.; Kluwer: Dordrecht, 1993; Chapter 13. (3) Rossky, P. J.; Simon, J. D. Nature (London) 1994, 370, 263. (4) Refs 5-8 contain reviews of various aspects of solvation dynamics. (5) Maroncelli, M. J. Mol. Liq. 1993, 57, 1. (6) Maroncelli, M.; Kumar, P. V.; Papazyan, A.; Horng, M. L.; Rosenthal, S. J.; Fleming, G. R. In Ultrafast Reaction Dynamics and SolVent Effects, Gauduel, Y., Rossky, P. J., Eds.; American Institute of Physics: New York, 1994; pp 310-333. (7) Bagchi, B.; Chandra, A. AdV. Chem. Phys. 1991, 80, 1. Bagchi, B. Annu. ReV. Phys. Chem. 1989, 40, 115. (8) Barbara, P. F.; Jarzeba, W. AdV. Photochem. 1990, 15, 1. Jarzeba, W.; Walker, G. C.; Johnson, A. E.; Barbara, P. F. Chem. Phys. 1991, 152, 57. (9) Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1988, 89, 5044. (10) Karim, O. A.; Haymet, A. D. J.; Banet, M. J.; Simon, J. D. J. Phys. Chem. 1988, 92, 3391. (11) Bader, J. S.; Chandler, D. Chem. Phys. Lett. 1989, 157, 501. (12) Belhadj, M.; Kitchen, D. B.; Blair, J. T.; Krogh-Jespersen, K.; Levy, R. M. J. Phys. Chem. 1991, 95, 1082. Levy, R. M.; Kitchen, D. B.; Blair,

Skaf and Ladanyi J. T.; Krogh-Jespersen, K. J. Phys. Chem. 1990, 94, 4470. (13) Rao, M.; Berne, B. J. J. Phys. Chem. 1981, 85, 1498. (14) Bursulaya, B. D.; Zichi, D. A.; Kim, H. J. J. Phys. Chem. 1995, 99, 10069; J. Phys. Chem. 1996, 100, 1392. (15) Muin˜o, P.; Callis, P. R. J. Chem. Phys. 1994, 100, 4093. (16) Ando, K.; Kato, S. J. Chem. Phys. 1991, 95, 5966. (17) Zhu, J.; Cukier, R. I. J. Chem. Phys. 1993, 98, 5679. (18) Fonseca, T.; Ladanyi, B. M. J. Phys. Chem. 1991, 95, 2116. (19) Fonseca, T.; Ladanyi, B. M. J. Mol. Liq. 1993, 60, 1. (20) Phelps, D. K.; Weaver, M. J.; Ladanyi, B. M. Chem. Phys. 1993, 176, 575. (21) Brown, R. J. Chem. Phys. 1995, 102, 9069. (22) Kumar, P. V.; Maroncelli, M. J. Chem. Phys. 1995, 103, 3038. (23) Skaf, M. S.; Ladanyi, B. M. J. Mol. Struct. 1995, 335, 181. (24) Maroncelli, M. J. Chem. Phys. 1991, 94, 2084. (25) Ladanyi, B. M.; Stratt, R. M. J. Phys. Chem. 1995, 99, 2502. (26) Carter, E. A.; Hynes, J. T. J. Chem. Phys. 1991, 94, 5961. (27) Ladanyi, B. M.; Skaf, M. S. Annu. ReV. Phys. Chem. 1993, 44, 335. (28) Olender, R.; Nitzan, A. J. Chem. Phys. 1995, 102, 7180. (29) Ladanyi, B. M.; Skaf, M. S. J. Phys. Chem. 1996, 100, 1368. (30) Ferrario, M.; Haughney, M.; McDonald, I. R.; Klein, M. L. J. Chem. Phys. 1990, 93, 5156. (31) Skaf, M. S.; Ladanyi, B. M. J. Chem. Phys. 1995, 102, 6542. (32) Tanaka, H.; Gubbins, K. E. J. Chem. Phys. 1992, 97, 2626. Koh, K. A.; Tanaka, H.; Walsh, J. M.; Gubbins, K. E.; Zollweg, J. A. Fluid Phase Equilib. 1993, 83, 51. (33) (a) Pa´linka´s, G.; Hawlicka, E.; Heinzinger, K. Chem. Phys. 1991, 158, 65. (b) Pa´linka´s, G.; Bako´, I.; Heinzinger, K.; Bopp, P. Mol. Phys. 1991, 73, 897. (34) Bertolini, D.; Cassettari, M.; Salvetti, G. J. Chem. Phys. 1983, 78, 365. Bertolini, D.; Cassettari, M.; Salvetti, G.; Tombari, E.; Veronesi, S. J. Non-Cryst. Solids 1991, 131-133, 1169. Mashimo, S.; Kuwabara, S.; Yagihara, S.; Higasi, K. J. Chem. Phys. 1989, 90, 3292. (35) Zhang, X.; Jonah, C. D. J. Phys. Chem. 1996, 100, 7042. Lin, Y.; Jonah, C. D. In Ultrafast Dynamics of Chemical Systems, Simon, J. D., Ed.; Kluwer: Dordrecht, 1994; pp 137-162. Lin, Y.; Jonah, C. D. J. Phys. Chem. 1993, 97, 295. Lin, Y.; Jonah, C. D. J. Phys. Chem. 1992, 96, 10119. (36) Decle´my, A.; Rullie`re, C.; Kottis, P. Laser Chem. 1990, 10, 413. Decle´my, A.; Rullie`re, C. In Ultrafast Reaction Dynamics and SolVent Effects; Gauduel, Y.; Rossky, P. J., Eds.; American Institute of Physics: New York, 1994; pp 275-295. (37) Yu, J.; Berg, M. Chem. Phys. Lett. 1993, 208, 315. (38) Chapman, C. F.; Fee, R. S.; Maroncelli, M. J. Phys. Chem. 1995, 99, 4811. (39) Carter, E. A.; Hynes, J. T. J. Phys. Chem. 1989, 93, 2184. (40) See, for example: Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (41) See, for example: Rappe´, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (42) Haughney, M.; Ferrario, M.; McDonald, I. R. J. Phys. Chem. 1987, 91, 4934. (43) Jorgensen, W. L.; Chandrasekhar, J.; Impey; R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (44) Fonseca, T.; Ladanyi, B. M. J. Chem. Phys. 1990, 93, 8148. (45) Skaf, M. S.; Fonseca, T.; Ladanyi, B. M. J. Chem. Phys. 1993, 98, 8929. (46) Neumann, M. J. Chem. Phys. 1986, 85, 1567. (47) Luzar, A.; Chandler, D. J. Chem. Phys. 1993, 98, 8160. (48) Hockney, R. W. Methods Comput. Phys. 1970, 9, 136. (49) Ciccotti, G.; Ryckaert, J.-P. Comput. Phys. Rep. 1986, 4, 345. (50) Ciccotti, G.; Ferrario, M.; Ryckaert, J.-P. Mol. Phys. 1982, 47, 1253. (51) Chandra, A.; Bagchi, B. J. Chem. Phys. 1991, 94, 8367.

JP961634O