2116
J. Phys. Chem. 1991, 95, 2116-2119
Breakdown of Linear Response for Solvation Dynamics in Methanol Teresa Fonseca* and Branka M. Ladanyis Department of Chemistry, Colorado State University, Fort Collins, Colorado 80523 (Received: December 1I , 1990)
We present the results of a molecular dynamia (MD) simulationstudy of solvation dynamics in methanol. Using nonequilibrium MD, we have calculated the relaxation rate of the solute-solvent potential energy change due to an instantaneous dipole creation in a diatomic solute, a model for solvent response to a charge-transfer electronic transition of a probe molecule. Linear response (LR) approximations to this process were obtained from equilibrium MD. LR results for methanol in the presence of the ground- and excited-state solutes differ substantially from each other and from nonequilibrium solvation dynamics. We show that 0-H bond motions dominate early solvation dynamics, but that other relaxation mechanisms subsequently become important, and explain why LR breaks down much more seriously for methanol than for other polar solvents previously studied by MD.
Introduction The dynamics of solvation has been the subject of increasing experimental'-' as well as theoreticalS13 study. Experiments measure the time-dependent Stokes shift of the fluorescence emission of probe molecules dissolved in polar solvents. The basic concept underlying the interpretation of the experimental results, in terms of the dynamics of the solvent molecules, is outlined in Figure 1. There, the dependence of the nonequilibrium free energies of the solvated ground and excited states of the probe is shown as a function of a generalized solvent ~ m r d i n a t e . ' ~ Optical excitation of the probe changes its charge distribution and produces, by the Franck-Condon principle, a nonequilibrium configuration of the solvent molecules around the probe. Subsequent relaxation of the solvent is accompanied by a time-dependent fluorescence spectral shift of the probe's emission; this shift is then used to quantify the solvent dynamics.14 Accompanying the intense experimental activity there have been a few molecular dynamics (MD) simulations of solvation dynamics. These studies have been performed in ~ a t e r , ' ~ a- ' ~ diatomic dipolar aprotic solvent,20and acetonitrile.21 Although the simulated systems are still far from the complex experimental ones, especially with respbct to the representation of the molecular probes, they are useful in helping us understand the fine microscopic details of the solvation process, and they constitute a critical test of the theories for solvation dynamics.S13 An interesting result, common to all the computational studies performed so far, is the fact that linear response (LR) theory applied to the equilibrium dynamics of the solvent in the presence of the ground state of the probe is able to predict remarkably well the major portion of the nonequilibrium response even when the probe undergoes a large change in its charge distribution.'6J8~20*21 In this Letter, we present a preliminary account of the results of a MD study of solvation dynamics in methanol. In contrast to previous studies, our results show that for methanol LR constitutes a very poor approximation to the nonequilibrium results. To our knowledge, this constitutes the first example of a clear breakdown of LR in solvation dynamics. Models and Methods We will consider two different model systems, each consisting of a solute molecule dissolved in methanol. Both model solutes are rigid diatomic molecules that differ only in size, Le., atomic diameters. The solute ground state is nonpolar (NP) and is represented by uncharged atomic sites. In the excited state an electron is transferred between the two sites, creating a dipolar (DP) solute. For both solutes, the two sites are identical and the molecular bond is equal to 1.43 A. Upon excitation the dipole moment of both solutes changes from 0 D (ground state) to 6.87 D (excited state). Authors to whom correspondence should be addressed.
0022-3654/91/2095-2116$02.50/0
TABLE I: Potential Parameters for the Solvent rad Solute sitec Me 0 H
(A) Solvent-Methanolb de (clkdlK 0.297 91.2 -0.728 81.9
0.43 1
0
u/A 3.861 3.083 0
(B) Diatomic Soluted
small large
0 0
*l
87.9
3.083
d21
87.9
4.200
g;"
"4 = site partial charge; c = site LJ ener = site LJ diameter. *Geometry: RMd = 1.43 A; hH = 0.945 LMeO-H = 108°53'. 'The H site is the hydroxyl hydrogen. each solute, the two atoms have the same masses and LJ potential parameters. Geometry: bond length = 1.40 A. Atom mass: 30mH.
Our model of methanol is a rigid three-site (CH3group, 0, and hydroxyl H) model whose structuralz2and dielectricf3 properties ( I ) Kosower, E. M.; Huppert, D. Annu. Rev. Phys. Chem. 1986,37, 127. (2) Barbara, P. F.; Jarzeba, W. Ado. Phorochem. 1990, 15, 1. (3) Simon, J. D. ACC.Chem. Res. 1988, 21, 128. (4) Maroncelli, M.; MacInnis, J.; Fleming, G. R. Science 1989,243,1674. (5) Mazurenko, Y.T.; Bakshiev, N. G. Opt. SpectraPC. 1970, 28, 490. (6) Bagchi, B.; Oxtoby, D.; Fleming, G. R. Chem. Phys. 1984,86,257. (7) van der Zwan, G.; Hynes, J. T. J. Phys. Chem. 19SS,89,4181. (8) Loring, R. F.; Mukamel, S.J . Chem. Phys. 1987,87, 1272. (9) Wolynw, P. G. J. Chem. Phys. 1987,86, 5133. (10) (a) Rip, I.; Klafter, J.; Jortner, J. J . Chem. Phys. 1988,88, 3246. (b) J . Chem. Phys. 1988,89,4288. (11) Nichols 111, A. L.; Calef, D. F. J. Chem. Phys. 1988, 89, 3783. (12) Zhou, Y.;Friedman, H. L.; Stell, G. J . Chem. Phys. 1989,91,4885. (13) (a) Bagchi, B.; Chandra, A. J. Chem. Phys. 1989, 90,7338. (b) Chandra, A,; Bagchi, B. J. Phys. Chem. 1989,93, 6996. (14) See, for example: ref 2 or Mataga, N. In Molecular In.teracrions; Vol. 2; Ratajczak, H., Orville-Thomas, W. J., Eds.; Wiley: New York, 1981. (15) Rao, M.; Berne, B. J. J. Phys. Chem. 1981.85, 1498. (16) Maroncelli, M.;Fleming, G. R. J . Chem. Phys. 1988, 89, 5044. (17) Karim, 0. A,; Haymet, A. D. J.; Banet, M.J.; Simon, J. D. J . Phys. Chem. 1988, 92, 3391. (18) Badcr, J. S.;Chandler, D. Chem. Phys. Lrrr. 1989, 157, 501. (19) Levy, R. M.; Kitchen, D. B.; Blair, J. T.; Krogh-Jesperscn, K. J. Phys. Chem. 1990, 94, 4470. Belhadj, M.; Kitchen, D. B.; Krogh-Jespersen, K.; Levy, R. M. J. Phys. Chem. 1991, 95, 1082. (20) (a) Hynes, J. T.; Carter, E. A.; Ciccotti, G.; Kim, H. J.; Zichi, D. A.; Ferrario, M.; Kapral, R. In Perspecrives in Phorosynrhesis; Jortner, J., Pullman, B., Eds.;Rcidcl: Dordrecht, 1990. (b) Carter, E. A.; Hynes. J. T. J . Chem. Phys., in press. (21) Maroncelli, M. J. Chem. Phys. 1991, 94, 2084. (22) Haughney, M.; Ferrario, M.; McDonald, I. R. J . Phys. Chem. 1987, 91, 4934. (23) Fonseca, T.; Ladanyi, B. M. J . Chem. Phys. 1990, 93, 8148.
0 1991 American Chemical Society
The Journal of Physical Chemistry, Vol. 95, No. 6, 1991 2117
Letters
where Ed is the electronic transition energy in the absence of the solvent and the overbar indicates an average in the presence of the perturbation that changes the solute electronic state at t = 0. In the MD simulation this average is taken over the different trajectories generated in the nonequilibrium runs. The experimental measurements are quantified in terms of the response function14 hv(t) - hv(m) C(t) = hv(0) - hv(m) and, using relation 2, C(r)can be calculated from the nonequilibrium MD results as C(t) = generalized solvent coordinate
Figure 1. Schematic illustration of the nonequilibrium free energy as a function of the generalized solvation coordinate for the solute ground (GS) and excited (ES) electronic states. The electronic excitation is shown as a full line with an upward arrow and the subsequent time-dependent fluorescence as dashed lines with downward arrows. u(0) and v(=) are the frequencies at the intensity maxima of absorption (and the initial, I = 0, fluorescence) and equilibrium fluorescence, respectively.
have been extensively studied. The intermolecular interactions are represented by sitesite Lennard-Jones (LJ) plus Coulombic potentials. A summary of the parameters used, for both the solvent and the solutes, is presented in Table I. The LJ interaction parameters between unlike sites have been evaluated by using the Lorentz-Berthelot combining rules. In the simulations we have used periodic cubic boundary conditions, with 255 molecules of methanol and one solute molecule in the original box whose side equals 25.9 A; the simulations were all performed in the NVE ensemble, with an average temperature of 298 K, and the equations of motion have been integrated by using SHAKE^) and leapfrog algorithms with a time step of 8 fs. The Ewald sum method has been used to treat the long-ranged electrostatic interaction^.^^ Two kinds of simulations have been performed for both solutes: first, equilibrium simulations where the solvent was equilibrated with the solute in either its ground or excited state. For these, 15 000-25 OOO time steps have been d e d in each case. Second, nonequilibrium simulations were performed in the following way: Initial configurations were extracted from the equilibrium runs for methanol in the presence of the ground-state solute. The solute state was then changed instantaneously from the ground to the excited state, and the system was then allowed to equilibrate to this new condition. About 110 trajectories, separated by 0.8 ps, were recorded in each case. In this Letter, we test the validity of the LR theory to predict the nonequilibrium results. A full analysis of the diffmnces in relaxation observed for the two solutes will be presented elsewhere.
Results and Discussion A convenient microscopic variable to gauge the solvent relaxation is AEl6tl8-21J6
AE
=:
V&W
- V&OS
(1)
where V, is the soluttsolvent interaction energy with the solute in the excited state (ES)and ground state (GS). As the Franck-Condon principle applies to the electronic transition, AE is also a free energy quantity, and its relation to the experimental measurements, specifically to the frequency Y at the maximum of the fluorescence intensity, becomes clear if we go back to Figure 1
hv(t) = E(?) + E,,
(2)
(24) Ciccotti, G.; Ryckaert, J.-P. Compur. Phys. Rep. 1986,4, 345. (25) Allen. M. P.; Tildesley, D. J. Compurer Slmulurion of Liquids; Clarendon Rcs: Oxford, 1987. (26) Warshell, A. J . Phys. Chem. 1982,86, 2218.
D(t)- G ( m ) m(0)- D(-)
Linear response can be applied to approximate the nonequilibrium response function C(t) by S(t), the equilibrium time autocorrelation function of the fluctuations in AE.27 S(t) can be calculated for the solvent in the presence of either the ground-state or the excited-state solute. In the former case, one obtains (4)
where (...)NP denotes an equilibrium ensemble average in the presence of the nonpolar, ground-state solute. C(t) is then a p proximated by use of the solvent dynamics in the absence of the solute dipole. Note, however, that the ensemble average in eq 4 uses the solvent configurations that constitute starting points for the nonequilibrium simulations. Alternatively, C(t) can be approximated in terms of the dynamics representative of the solvent in equilibrium with the excited-state, dipolar solutez0
where ( . . . ) D p denotes an equilibrium ensemble average in the presence of the excited probe and 6AE = AE - (AE)Dp.(The average value of AE is nonzero when the solute is in the excited state.) Since the solut-olvent interactions are uite different for the ground- and excited-state solutes, SoS(?),S (t), and C(t) might be expected to differ from each other. However, until now, all the MD studies of solvation dynamics have shown that Sos(?)was able to predict most of the time evolution of C(t) remarkably we11,16J8.m.21 even though the changes involved in the excitation of the probe were rather large and, therefore, poor agreement was expected. Using our simulation results for methanol, we have evaluated C(t) and its equilibrium analogues, Sos(?)and Ses(t), in the presence of both the small and large solutes. The comparison between the LR results and C(r) is presented in Figure 2. As can be seen from these results, LR works quite poorly in each case although the agreement is somewhat better in the case of the larger solute. The last is consistent with fact that the perturbation inflicted on the solvent is smaller with the large solute. Despite the poor agreement, several interesting facts emerge from the analysis of our results, and we discuss them in what follows. For both cases, the initial Gaussian decay of C ( t ) - o f about 20 fs duration-is well predicted by P ( t ) . (Aswe discuss below, the Gaussian character of the initial decay of C(r) is, in fact, responsible for the good agreement between LR and the nonequilibrium solvation results in other solvents.28) In view of its
Bs
(27) We u8c the symbol c(r) for the nonequilibriumresponse function in order to make our notation consistent with that of most of the experimental papers on timedependcnt fluorescence. However, readers should note that other symbols have been used for this function in recent MD simulation work, in order to distinguish it from the yilibrium time correlation, here denoted by S(r). In some of this work,16.202the symbol S(r) was used to denote the nonequilibrium response function. (28) The importance of the initial inertial dynamics of c(r)was first pointed out in ref 20.
2118 The Journal of Physical Chemistry, Vol. 95, No. 6, 1991
Letters
(a) Small solute
, - O . l 0.8
1's.
I
*
0.0
l
0.2
_____
,-
~
I
0.4
0.6
time(ps)
-_
~
I
~
-
~
(b) large solute
1.0
'I:, '\ ' ,
I
20
sGyri
40
60
80
time(fs)
Figure 3. Effects of deuteration of the hydroxyl hydrogen of methanol on the time evolution of the nonquilibrium res nse function C(r) following dipole creation in the large ( u = 4.200 cdiatomic solute. The results for CH30H solvent are depicted by a full line and those for the CH30D solvent by a dashed line. 1-0.1
t\ o*9t
i
\
" t \ V
. -.-._.--_._____. -..
0.0
0.0
0
I
0.8
0.2
0.4
0.6
0.8
time(ps)
-6.0
Figure 2. Comparison of the&ne evolution of the non-&uilibrium response function, C(t) (full line), to the linear response time correlations for methanol in the presence of the ground-state, nonpolar solute, P ( r ) (dash-dotted line) and in the excited-state, lar solute, Ses(t)(dashed line). The results for the small ( u = 3.083 and large ( u = 4.200 A) diatomic solutes are depicted respectively in parts a and b of the figure.
5
importance, we have attempted to identify the types of motions responsible for this initial decay. Considering the time scale involved-a few femtodecades-our educated guess would suggest motions that involve the hydrogen atom of the hydroxyl group in methanol: these motions can be either the free rotation of the 0-H bond around the M e 4 bond (the corresponding moment of inertia is about 10 times smaller than those involving rotation of the heavy particles) or the hindered librational motion along the intermolecular hydrogen bond that in effect corresponds to frustration of the rotation just described. To investigate this assumption and to distinguish between the contributions of free rotation and librational motion to the nonequilibrium response, we have performed several special nonequilibrium simulations that we now describe. In an attempt to identify the time scale associated with librational motions, we have redone our nonequilibrium simulations in deuterated methanol, CHpOD. A comparison between the nondeuterated and deuterated results for the large solute is shown in Figure 3. (Similar results were obtained for the small solute.) From Figure 3 we can identify the shoulder in C(t) a t about 40 fs in nondeuterated methanol as due to librational motion-the frequency of this bounded motion is shifted by 4 2 toward lower values in CH30D. Note the important fact that this time scale is late in the initial decay. To investigate the role of free 0-Hbond rotation, we performed nonequilibrium simulations in which we increased 10-fold the oxygen and methyl group masses of the solvent molecules that constitute the first solvation shell of the solute in its ground state: in this way, all motions but the free H rotation are slowed down considerably; indeed, they are practically frozen during the simulation time. Figure 4 shows the comparison between C(t) for 'heavy methanol" and the normal one in the presence of the large solute. As Figure 4 clearly shows, the Gaussian decay of C(t)
0
20
40
60
80
time(ft)
Figure 4. Time evolution of the nonquilibrium response function C(r) following dipole creation in the large (u = 4.200 A) diatomic solute. The contribution of the 0-H bond dynamics to early relaxation is illustrated by comparing the results for ordinary (full line) and "heavy" (dashed line) methanol. In the latter solvent the masses of the Me and 0 sites have been increased 10-fold.
remains unchanged. This is an extremely important result, as it shows that the initial decay of C(t) is determined by the free H rotation about the M A bond. (Moreover, this is also responsible for the different dynamics exhibited by methanol in the presence of the two solutes. As shown in a forthcoming work, the large solute has a structure-breaking effect on methanol, providing an increased number of freely rotating 0-H bonds that relax C(t) more effectively than for the small solute.) We have thus identified the inertial motion responsible for the initial decay of C(t),and we can now discuss the origin of the breakdown of LR for methanol, as well as its success for other solvents, in particular for water. To start our diaussion, let us rephrase the ideas behind LR: if the perturbation imposed on the system is small enough, then we can describe the response of the system using its statics and dynamics in the absence of the perturbation. However, the perturbations we have been discussing so far are not small at all: there is a net change of 6.87 D in the solute dipole moment. Evidence that this perturbation is indeed viewed by the solvent as rather large is provided by the completely different equilibrium dynamics of methanol in the presence of the GS and Es of the probe (see Figure 2). Different equilibrium dynamics were also observed for the other solute-solvent systems,16-20.2'Therefore, the only thing that makes sense for us to ask is, which part of the solvent dynamics is not strongly affected by the perturbation? Just after the perturbation, different motions of the solvent molecules will be animated with the velocities they had just before the perturbation; inertia will make them preserve these motions until interactions, e.g., restoring forces from hydrogen bonds or short-range repulsions, begin to take effect. Thii inertial motion is not affected by the perturbation, and LR should
J. Phys. Chem. 1991,95. 21 19-2121 correctly describe this part of the relaxation (even for large perturbations), if the average weights are taken from the appropriate equilibrium configurations, i.e., those used as starting points for the nonequilibrium response. (See also ref 20.) Now, the question we should ask is, how much of C(t) is due to this inertial motion? According to our discussion, it is clear that LR will work if most of C(t)is determined by inertial motion: indeed, LR should accurately describe this part of C(t). For the solvents studied so far we have the following: (i) For the diatomic dipolar solventz0and the three-site model of acetonitrile2' (no hydrogen bonding), the initial Gaussian decay of C(t) constitutes 7O-80% of the total relaxation and occurs on a time scale of 100-250 fs. SGS(t)accurately describes this part of C(t).20,2'This time scale is considerably longer than that characteristic of hydroxylic solvents (see below) because there are no small moments of inertia here.29 The effectiveness of inertial dynamics in relaxing C(t) can be ascribed to the fact that there are no specific interactions, and librational dynamics is thus unimportant. (ii) For methanol and water,l6,l8 pS(t) still describes the Gaussian decay of C(t),which now occurs on a much shorter time scale (-30 fs) than in aprotic solvents. Librations, resulting from the hydrogen-bond restoring forces, come into play at an early stage of the relaxation and break this initial inertial decay. However, we still have to understand why water is so different from methanol: in water we observe that the Gaussian decay of C(r)constitutes 50-80% of the total response'6J8 while in methanol it accounts for only 1O-20% (see Figure 2). We have shown here (29) In the model of acetonitrile used in ref 21, the methyl group is represented by a single interaction site, so the molecule is effectively linear.
2119
that the initial decay of C(t) for methanol is determined by free H rotation, in other words the nonfrustrated part of librational motion. Although such a precise identification has not been made for water, we believe that the same is true there. The key difference is that in water there are three distinct ways of performing this rotational motion (with nearly equal moments of inertia), while in methanol there is only a single way. Therefore, the effectiveness of the free H rotation in relaxing C(t)should be much larger in water than in methanol-this is just what is observed. Finally we examine the long time behavior of C(t). We see that it is quite similar to the long time dynamics predicted by the LR of the excited state, SES(t)(see Figure 2). We also note that the long time relaxation rate of C(t) is virtually independent of solute size: this is an indication that these dynamics are determined by solvent molecules far away from the solute and for which the perturbation is effectively small as it is shielded by the molecules in the first solvation shell. We can thus understand why SES(t)predicts well the long time decay of C(t). This fact has first been pointed out by Carter and Hynes for their diatomic solventzo and is also observed in acetonitrile2' and water.I6 Summarizing, we have observed and discussed here the limitations of LR in describing solvation dynamics in methanol. For the first time, we have identified the motions responsible for the initial Gaussian decay of C(t)in methanol and presented the first rationale for the differences in solvation dynamics of methanol and water.
Acknowledgment. This work was supported in part by N I H Grant GM27945 and NSF Grant CHE-8808206. Acknowledgment is also made to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research.
A *H NMR Study of the Formation of Ethylldyne from Acetylene and Hydrogen Coadsorbed on Pt C.A. Klug? C . P . Slicl~ter,*~~~*~s Departments of Physics and Chemistry and Materials Research Laboratory, I 1 IO W. Green Street, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801
and J. H.Sinfelt Exxon Research and Engineering Company, Annandale, New Jersey 08801 (Received: January 7, 1991)
The authors report the use of deuterium NMR to study the conversion of vinylidene (CCH2) plus hydrogen to ethylidyne (CCHp)using small supported Pt metal clusters typical of industrial catalysts. They find an activation energy of 9.3 f 2.5 kcal/mol and a preexponential factor of 3.8 X s-*. They show that the reaction explains previous observations of some carbon-carbon bond scission in the temperature region of 450-650 K for acetylene adsorbed on Pt.
Introduction Some of the earliest studies of molecules adsorbed on metal surfaces involved determining the room-temperature structures of the simple hydrocarbons acetylene and ethylene adsorbed on Pt. It is generally agreed that near room temperature ethylene loses a hydrogen atom and converts to ethylidyne (CCHp)while acetylene rearranges to form vinylidene (CCH2). Several researchers including Somorjai' and Ibach2 found that background hydrogen in their ultrahigh-vacuum experimental apparatus had allowed the vinylidene to convert to ethylidyne, which led them *To whom correspondence should be addressed. 'Department of Physics. f Department of Chemistry. t Materials Research Laboratory. 0022-3654/91/2095-2119%02.50/0
to propose initially that both acetylene and ethylene formed ethylidyne. We were led to our present work by earlier studies in our group of carbon-carbon bond scission for CCH3 and CCH2 For CCH3 adsorbed on either Pt or Ir, Wang et aL3 found that the carbon-carbon bond scission took place at 450 K. For CCH2, however, the results were more complicated. Whereas, for Ir, we found that after acetylene was adsorbed at room temperature the carbon-arbon bonds broke at a temperature of about 450 K, for Pt, we found a substantially higher temperature of carbon-carbon (1) Kesmodel, L. L.; Dubois, L. H.; Somorjai, G. A. Chem. Phys. Leu. 1978, 56, 261. (2) Ibach, H.; Hopster, H.; Sexton, B. Appl. SurJ Sci. 1977, 1, 1. (3) Wang, P.; Slichter, C. P.; Sinfelt, J. H. J. Phys. Chem. 1990, 94, 1154.
0 1991 American Chemical Society
