Solvation Dynamics of C153 in Supercritical Fluoroform: A Simulation

May 2, 2006 - Molecular dynamics (MD) simulations of a probe solute (coumarin C153) in supercritical fluoroform are used to study time-dependent ...
0 downloads 0 Views 228KB Size
Download PDF
10120

J. Phys. Chem. B 2006, 110, 10120-10129

Solvation Dynamics of C153 in Supercritical Fluoroform: A Simulation Study Based on Two-Site and Five-Site Models of the Solvent Francesca Ingrosso* and Branka M. Ladanyi Department of Chemistry, Colorado State UniVersity, Fort Collins, Colorado 80523 ReceiVed: February 23, 2006; In Final Form: March 30, 2006

Molecular dynamics (MD) simulations of a probe solute (coumarin C153) in supercritical fluoroform are used to study time-dependent solute-solvent interactions. We study the dynamics of solvent reorganization in response to electronic excitation of C153 at a temperature of 1.03 Tc (the critical temperature) and a series of densities above and below the critical density. Simulations of a two-site and five-site models of fluoroform are presented and compared. The time-dependent solvation response after solute electronic excitation is studied in the two cases, and the five-site results present an earlier onset of exponential decay that is closer to what is expected to be the experimental response. This is confirmed by comparison to experiment. In addition to obtaining the solvation response from nonequilibrium MD trajectories, approximate solvation responses were obtained from equilibrium time correlations of the fluctuations in the solvation energy change in the presence of ground- and excited-state solutes. For the five-site model, the equilibrium excited-state response shows stronger density dependence than the ground-state one. The nonequilibrium response appears to have an intermediate decay rate between the two equilibrium functions. The solute-partial-charge-solvent-induceddipole interaction was also taken into account by means of a perturbative approach, which improved the agreement with experimental measurements available at densities corresponding to 1.4-1.6 Fc (where Fc the critical density). From the comparison between the two models, it is possible to conclude that an atomistic description is necessary for correctly representing the portion of solvation dynamics that is related to reorientation. This consideration is supported by providing results for orientational time correlation functions and by comparing the correlation times with the experimental ones.

1. Introduction Fluids in the supercritical (SC) regime experience intermediate behavior between gas-phase-like (low-density limit) and liquidphase-like (high-density limit). The SC regime is characterized by density inhomogeneities, and small variations of the system pressure produce changes in the bulk density due to the high compressibility in the near-critical region.1,2 This property can be exploited to tune the solvent effects on reactions taking place in the SC regime.1,3 Because of the interest in chemical reactivity, most studies in the supercritical regime are now focusing on augmentation effects related to the presence of an “attractive” solute (i.e., a solute whose interaction with the solvent is stronger than the interaction between solvent molecules) and their relation with solute-solvent interactions.4-8 Our interest here is focused on a particular aspect of solvation in SC fluoroform, namely the dynamics of the solvation process. To avoid large density fluctuations close to the critical point, we performed our study in the near-critical regime, where macroscopic scale fluctuations are not present but local density augmentation is observed, and from now on the discussion will concern this portion of the phase diagram. Solvation dynamics in SC fluids was investigated by Egorov and collaborators6-10 using theoretical and computer simulation methods. A nonlinear generalized Smoluchowski-Vlasov equation and mode-coupling theory were used to describe the * Corresponding authors. E-mail: [email protected] (B.M.L.).

[email protected]

(F.I.);

dynamical aspects of solvation in SC fluids, with a particular interest in studying the differences between nonequilibrium and equilibrium results. Simulations in SC carbon dioxide11 and water12 performed with an atomistic model for the solvent were carried out by other groups. Results from theory and simulations showed that the slow portion of the decay of the solvation response function is influenced by the changes in the local density, and that breakdown of linear response theory is observed for moderate and low densities in the SC regime. This is to be expected because deviations from bulk behavior are associated with density inhomogeneities, and these deviations become smaller as the density increases toward liquidlike values. Theoretical/computational investigations of time-dependent solvation response in SC fluids have been carried out on a variety of different systems (see refs 6-15 and ref 1 for a review of older literature). The work reported here is based on an atomistic representation of the solute and both on a two-site and on an atomistic description of the solvent, and we show here that the latter is necessary to achieve an accurate description of the librational-diffusive portion of solvation dynamics. To conduct this study, we performed molecular dynamics (MD) simulations of coumarin C153 (see Figure 1) in two-site16 and five-site17 fluoroform. C153 is one of the most common polar solvation probes in experiments18 based on the consideration19 that the S0 f S1 transition is associated with a large increase in the molecular dipole moment, with a measurable effect on polar solvents. As for the solvent, fluoroform is a dipolar SC fluid that interacts more strongly with the solute than do other room-temperature SC fluids such as CO2 and

10.1021/jp061170m CCC: $33.50 © 2006 American Chemical Society Published on Web 05/02/2006

Solvation Dynamics of C153 in Supercritical Fluoroform

J. Phys. Chem. B, Vol. 110, No. 20, 2006 10121 shift is due to the change, ∆E(t), in the solute-solvent potential23 that occurs as a result of perturbing the preexisting solute-solvent equilibrium. The perturbation, ∆E, is turned on at t ) 0, and t ) ∞ represents equilibrium in the presence of the solute excited state. It is often assumed that ∆E is due solely to the change ∆qR in the solute partial charges.24 In this case, ∆E is electrostatic and can be expressed as:25,26

Figure 1. The solute, C153. Labeling used in the text for indicating atomic sites is shown.

ethane. This gives rise to larger values of the Stokes shifts, which make experimental measurements in this solvent more accessible than in other room-temperature SC solvents. Experimental measurements15 of solvation dynamics have been recently performed on a similar system (4-dimethylamino4-cyanostilbene (DCS) in SC fluoroform at densities higher than the critical). A comparison of our simulated result with experiment is therefore possible. By using pulsed laser experiments,20,21 it is possible to monitor the time evolution of the shift in the fluorescence spectrum. The function of interest, the time-dependent Stokes shift S(t), is defined as:

S(t) )

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

(1)

where ν(0) is the value of the maximum fluorescence band immediately after electronic transition, and ν(∞) is the corresponding value at equilibrium. S(t) gives information on the solvent response induced by a sudden change in the solute charge distribution. In Section 2, we define the theoretical counterpart of this function, as we calculated it from MD trajectories. Calculations of S(t) were performed for the two-site and the five-site models. We ran equilibrium simulations with the solute in the ground and in the excited state, and a smaller group of densities was chosen to run nonequilibrium simulations. Differences in equilibrium and nonequilibrium results were analyzed, as well as differences between ground- and excited-state equilibrium responses. In addition, we are interested in investigating the role of induction in solvation dynamics. We already found5 that induction makes a nontrivial contribution to effective local density in the near-critical regime. In the present work, we calculated the solvation response of fluoroform at 1.5 Fc by including the contribution due to interaction between the solute charge distribution and induced dipoles on the solvent molecules to the interaction potential by means of first-order perturbation theory. This density was chosen in order to compare with experiments,15 which were run at densities larger than 1.4 Fc. Finally, we extended our study to orientational dynamics in pure fluoroform. We calculated orientational time correlation functions for the five-site model and compared the rotational correlation times to experimental measurements available in the literature.22 A comparison between single-molecule and collective orientational dynamics in the two-site and in the five-site models was also performed. This paper is organized as follows. In Section 2, the theoretical concepts used in this study are presented, followed by a description of the computational procedures adopted in Section 3. Results are reported and discussed in Section 4, while we summarize the major findings of this work in Section 5. 2. Theoretical Background Under the assumption that the variation of the solute geometry in the excitation process is negligible, the fluorescence Stokes

∑

∆E )

∆qRVR,

(2)

R∈solute

The electrostatic potential VR at the solute site R is defined as:

V(rR,jβ) )

Nmol

1 4π0

qjβ

∑ ∑ j)1 β∈solvent r

(3)

R, jβ

where qjβ is the partial charge on site β of the jth solvent molecule, rR, jβ is the distance between the solute site R and site β on the jth solvent molecule, and Nmol is the number of solvent molecules. The theoretical counterpart of S(t), defined in eq 1, is the nonequilibrium solvation response function (SRF), defined as:24

SRFneq(t) )

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

(4)

where the overbar denotes a nonequilibrium average in the perturbed system. This definition, based on nonequilibrium simulations, can be very demanding from a computational point of view for large systems. Ground-state equilibrium trajectories need to be saved at long enough intervals to ensure that starting configurations are independent. The saved configurations are then used as starting configurations for the nonequilibrium runs, and the results are averaged over the total number of the (nonequilibrium) trajectories. A large number of nonequilibrium trajectories is generally necessary for good statistics. If a perturbation of a system from equilibrium is small, its effects can be approximated in terms of the linear response theory, which connects the properties of the perturbed system to fluctuations in the unperturbed system.27 In the present context, the response of the solution to the perturbation represented by ∆E(t) can be approximated in terms of the time correlation of δ∆E ) ∆E - 〈∆E〉 in the equilibrated unperturbed system (i.e., the solvent in the presence of the ground (0) or excited (1) state solute):24

SRFeq0,1(t) )

〈δ∆E(0) δ∆E(t)〉0,1 〈|δ∆E|2〉0,1

(5)

where 〈‚‚‚〉0,1 indicates that the equilibrium ensemble average is performed either in the ground or in the excited state. The range of validity of linear response for solvation dynamics can be tested by comparing SRF0eq(t) and SRF1eq(t) and both equilibrium responses to SRFneq(t). To investigate the effects of induction on solvation dynamics, we introduced the interaction between the solute charge distribution and the solvent-induced dipoles by using an approach based on first-order perturbation theory.5 This interaction gives rise to a term to be added to eq 2 when calculating the solute-solvent interaction energy ∆E. We used the equilibrium response (ground and excited state). The resulting ∆E was used in the calculation of SRF0,1 eq (t) at 1.5 Fc.

10122 J. Phys. Chem. B, Vol. 110, No. 20, 2006

Ingrosso and Ladanyi

We described the effect of the polarization on each solvent molecule j as due to the presence of an electric field ej generated by the solute, which gives rise to an induced dipole µind, j. The first-order perturbation to the change in the total interaction energy can be written as Nsolv

∆Eind ) -

1 1 0 0 [µind, ∑ j·ej - µind, j·ej ] j)1

(6)

where enj is the electric field generated by the solute charge distribution at the center of mass of the solvent molecule j. The solvent-induced dipole for the nth solute can be expressed in n n terms of the solvent polarizability tensor rj: µind, j ) rj‚ej (further details are contained in ref 5). The solvation response function describes the decay of the solvent polarization as a consequence of a change in the solute charge distribution. In the short time scale, the behavior of this function is related to inertial/reorientational dynamics in the fluid.20 Solvent reorientational dynamics can be described in terms of time correlation functions of vectors embedded in the molecules. Among these vectors, the molecular dipole is the most relevant to the electrostatic solvation response. We consider the single-molecule orientational time correlation functions:

Ψl(t) ) 〈Pl[uˆ (0)·uˆ (t)]〉

(7)

where uˆ is a unit vector along the direction of the molecular dipole, and Pl is the Legendre polynomial of rank l. The correlation times τl associated to these functions are defined as:

τl )

∫0∞Ψl(t) dt

(8)

The first rank collective dipole time correlation function N µi (where µi is the (TCF) is defined in terms of M ) ∑i)1 single-molecule dipole moment, and N is the total number of molecules):

Ψ1c(t) ) 〈M(t)·M(0)〉

(9)

and the corresponding correlation time is obtained from the normalized TCF as:

τ1c )

Ψ (t)

∫0∞ Ψ 1c(0) dt

(10)

1c

3. Computational Details Simulations of all-atom coumarin C153 in two-site fluoroform16 and in all-atoms fluoroform17 were performed by using a locally modified version of DL POLY 2.28 Four reduced densities (F/Fc ) 0.13, 0.59, 0.95, and 2.0) were considered in the two-site case, and five reduced densities (F/Fc ) 0.13, 0.3, 0.59, 1.5, and 2.0) were considered in the five-site case. Simulations were run in the microcanonical ensemble at a temperature slightly higher than the critical (T/Tc ) 1.03) to avoid macroscopic density fluctuations. We briefly report some features of the two models. Further details are collected in the original papers, refs 16 and 17. In Table 1 experimental and calculated critical points are collected. The two-site model16 was designed to accurately reproduce thermodynamic and static dielectric properties of fluoroform in the SC regime with an affordable computational effort. The two sites are identical except for the charge, and the masses (m

TABLE 1: Critical Variables from Experiment and from Simulation experimental43 five-site17 two-site16

Fc (g/cm3)

Tc (K)

0.526 0.630 0.541

299.1 309.0 300.2

) 35.01 amu) and separation (d ) 1.670 Å) were chosen to reproduce the rotational constant of real fluoroform. All the parameters of the model were tuned to reproduce experimental measurements (liquid-vapor coexistence, second pressure virial coefficient, viscosity, dielectric constant). The five-site model17 was developed to calculate orthobaric curves of fluoroform due to the interest for technological applications of this fluid as a refrigerant. The geometrical parameters were provided by accurate quantum calculations, and the potential was optimized to reproduce the liquid-vapor coexistence curve between 180 and 270 K and the pair correlation functions for the liquid from coherent neutron scattering experiments. The solute ground- and excited-state geometry and partial charges were obtained from quantum mechanical calculations based on the density functional theory. More details are contained in ref 5. In the simulation runs, the box was cubic and it included 255 solvent molecules and one solute molecule for densities lower than 0.95 Fc, while the system size was increased to 499 solvent molecules and one solute molecule at higher densities. The dimensions of the box were adapted to reproduce the experimental fluid densities16 at the temperature of interest, taking into account the larger size of the solute molecule compared to the solvent. We performed equilibrium simulations of ground- and excited-state C153 for each density and for both solvent models. In the two-site case, 500-800 independent trajectories were also saved from the ground-state equilibrium run at four different densities and used to perform nonequilibrium simulations, as described is Section 2. The same approach was used to calculate the nonequilibrium response for the five-site model at 0.3 Fc and 1.5 Fc. Averages were performed on a set of 800 trajectories. The time step used in all simulations was 2 fs. For both solvent models, equilibration required different times at different densities: longer equilibration time (4 ns) was necessary for densities equal to and smaller than 0.95 Fc, while 2 ns equilibration was run for higher densities. Data acquisition in the equilibrium simulations was performed on 4 ns trajectories at all densities but the lowest, for which we employed 5 ns trajectories. Nonequilibrium simulations were run for 4 ps. Simulations similar to those described for the solutions were run for the pure solvent, with a simulation box containing 256 and 500 fluoroform molecules below and above the critical density, respectively. The intermolecular potential was the combination of Coulomb and Lennard-Jones pair interactions between sites m and n, separated by a distance r, on different molecules:

umn(r) ) 4mn

[( ) ( ) ] σmn r

12

-

σmn r

6

qmqn + 4π0r

(11)

where the parameters σ and  are associated to the well depth and the diameter of the LJ potential, respectively. LorentzBerthelot combination rules were employed,29 with the exception of the LJ parameters for the interaction between solvent sites H and F (five-site model), which were tuned to better reproduce experimental data17 (see Table 2). The LJ potential at each

Solvation Dynamics of C153 in Supercritical Fluoroform

J. Phys. Chem. B, Vol. 110, No. 20, 2006 10123 TABLE 4: Stokes Shift Calculated from Simulationa

TABLE 2: Intermolecular Potential Parameters for Fluoroform σ (Å)

/κB (K)

q (e)

85.48

(0.275

10.0 54.6 40.0 22.36

0.106 0.440 -0.182

two-site “CHF3”

3.50 five-site

H C F HF

2.17 3.15 2.975 2.3

a

Fc

Stokes shift3 (cm-1)

0.13 0.30 0.59 1.5 2.0

800 1400 1600 1600 1850

The error associated with the calculation is 200 cm-1.

TABLE 3: Intramolecular Potential for Five-Site Fluoroform, Eq 12 F-C-F H-C-F

θ (deg)

κ (kJ mol-1 deg-2)

108.63 108.64

0.1457 0.0984

density was cut off at half the simulation box. The Coulombic and LJ parameters for the solvent were the same as in ref 16 for the two-site model and as in ref 17 for the five-site model. The parameters used for the intermolecular potential in eq 11 in the two cases are reported in Table 2, where we refer to the two-site model particles as “CHF3”. For the solute, we used the Lennard-Jones parameters from ref 19 and, as already mentioned, partial charges from ref 5. The interaction potential included some intramolecular flexibility. In the case of five-site fluoroform, bending is taken into account by using a harmonic potential:

κ U(θ) ) (θ - θ0)2 2

(12)

The parameters in eq 12 are reported in Table 3. The bond lengths were kept fixed. Stretching, bending, and torsional degrees of freedom in the nonaromatic portion of the solute molecule were included in the intramolecular potential (for details, see ref 26). As for the density at which we included the effects of induced-dipole interaction (F ) 1.5 Fc), as described in Section 2, we also ran ground- and excited-state equilibrium simulations for C153 in five-site fluoroform in which the solvent geometry was kept rigid (the intramolecular distances and angles were fixed to the equilibrium values reported in Table 3). In all simulations, the SHAKE30 algorithm was used to implement the geometrical constraints for the nonflexible parts of the solute and of the solvent intramolecular structures. The integration of the equations of motion was performed according to the Verlet leapfrog algorithm, and the long-range interactions were calculated using the Ewald method with conducting boundary conditions.29 4. Results 4.1 Solvation Dynamics. In this section, we compare the predictions for solvation dynamics obtained using two-site and five-site fluoroform models. Before we discuss time-dependent properties, we present a comparison between calculated and experimental4 Stokes shifts for C153 in SC fluoroform. Our results are collected in Table 4. The experimental result for the Stokes shift is approximately constant in the range 0.3-2.0 Fc, and the value of the constant is 1550 ( 100 cm-1. Within the errors associated with the measurements and with the simulation, our results are in reasonable agreement with experiments.

Figure 2. Solvation dynamics in two-site supercritical fluoroform in the density interval considered. Equilibrium and nonequilibrium results are shown.

We now discuss the results obtained with the two-site model for the solvation response at different densities, shown in Figure 2. In the first panel, the nonequilibrium responses are shown. The most evident features are the presence of a dip below zero at short t and a weak density dependence of the results. Very similar responses are obtained for the equilibrium ground- and excited-state simulations, shown in the middle and bottom panels, apart from a much slower decay for the lowest density in the excited-state result. Results from the five-site model are shown in the top panels of Figures 3 and 4. The bottom panels display fits of the calculated functions, which will be discussed later. The densities displayed are different from the ones reported in Figure 2, but the density range considered is the same. First of all, we notice that the pronounced dip at short times found for the two-site SRFs is not present in the five-site case, and for the latter, all functions show a more regular time evolution in the short and intermediate time scales. The relaxation mechanism in these time scales is related to reorientational motion of the solvent.31 The high-frequency oscillations in the calculated five-site responses, more evident in the excited case results, are due to solvent intramolecular bending vibrations, which were included in the interaction potential, as discussed in Section 3. This feature is not present in the results shown in Figure 2 because they were obtained with a rigid, two-site model. With the exception of the lowest-density result, the two-site SRFs exhibit small deviations from linear response behavior. We shall show later that the results obtained from the five-site

10124 J. Phys. Chem. B, Vol. 110, No. 20, 2006

Figure 3. Solvation dynamics in five-site supercritical fluoroform: equilibrium ground-state responses (top, full lines) and fits to eq 15 (bottom, dashed lines). The top panel inset shows the details of he short-time dynamics.

Ingrosso and Ladanyi especially at distances corresponding to the first solvation shell. Differences become smaller as the density increases. The position of the first peak is shifted toward smaller distances at the lowest density, and the peaks are broader, because this is a gas-phase-like density. As a consequence of larger effective local density, the peak height is larger at lower densities and decreases as the bulk density increases. Because of different clustering in the ground and in the excited states, the system under study is likely to give rise to deviations from the linear response approximation at low and intermediate densities. This consideration agrees with the fivesite results for the dynamic solvation response. We also report here solute-solvent orientational correlations at four different densities obtained with the five-site model. We computed the projection of the pair correlation function along one of the coefficients of the rotational invariant expansion:33

hmnl(r12) )

∫ h(12)Φmnl(12) dΩ1 dΩ2 ∫ [Φmnl(12)]2 dΩ1 dΩ2

(13)

where h(12) ) g(12) - 1 is the total pair correlation between molecules 1 and 2, Φmnl(12) is a rotational invariant function, while Ω1 and Ω2 are solid angles describing molecular orientations. In particular, we considered the correlation function h,110 which is related to the Φ110 coefficient, giving information about relative orientation of a pair of molecules

Φ110(12) ) uˆ 1‚uˆ 2

Figure 4. Solvation dynamics in five-site supercritical fluoroform: equilibrium excited-state responses (top, full lines) and fits to eq 15 (bottom, dashed lines). The top panel inset shows the details of the short-time dynamics.

model at two different densities exhibit larger deviations from linear response than those observed for the two-site model. In ref 11, the authors found that the nonequilibrium response shows a very similar behavior to the excited-state (equilibrium) one, whereas the ground-state response is slower at low and intermediate densities in the near-critical region. In addition, other simulation-based works8,12 predict a behavior that is more similar to what we found by employing the five-site model. The validity of the linear response approximation is related to differences in the solvation environment experienced by the solute in the ground and in the excited states.32 Local enhancement of solvent density around the solute occurs through translational diffusion to the first solvation shell. If clustering in the ground state differs from clustering in the excited state, larger deviations from linearity are observed when effective local densities are higher (low and intermediate densities in the nearcritical regime).11,8,12 To illustrate this issue, we present the results that we obtained for solute-solvent site-center-of-mass pair correlation functions for five-site fluoroform. We report in Figure 5 results for three densities (0.13, 0.3, and 1.5 Fc) and for six C153 sites (see Figure 1 for labeling). The plots show the details around the first peak. At the lowest density, the difference in solvent clustering around the ground- and the excited-state solute is large,

(14)

We defined the three unit vectors uˆ 1, uˆ 2, and rˆ12, in such a way that uˆ 1 is oriented along the dipole moment of a solvent molecule, uˆ 2 is oriented along an internuclear vector in the solute molecule, and rˆ12 is along the distance between the solvent center of mass and the center of mass of the solute sites considered. Within the solute, we chose two different uˆ 2 vectors, one along the C3 f C8 direction, the other one along the N13 f C2 direction (see Figure 1), and we analyzed the functions h110 for each of them. The reason for this choice, discussed in detail in ref 26, is related to the fact that the charge transfer from the donor (amino group) to the acceptor (CF3) group in C153 during excitation is well represented by considering these two internuclear vectors. Results are shown in Figure 6 for the C3-C8 direction and in Figure 7 for the N13-C2 direction. First, we observe the density dependence of the calculated functions along both internuclear directions. The first peak, occurring at distances corresponding to the first solvation shell, becomes less positive (ground state) or less negative (excited state) as the density increases. Second, there is a prominent difference in the solvation shell structure between the ground and the excited states. This difference is present at all densities, but more evident at the low and intermediate densities, where augmentation effects are larger. In the ground state, orientation of solvent molecules with respect to both solute internuclear vectors is preferentially parallel, on average, at distances corresponding to the first solvation shell, On the other hand, in the excited state, the antiparallel orientation is preferred. This effect is related to different polarity along the two C153 vectors in the ground and in the excited states.5 The five-site, equilibrium ground- and excited-state solvation response functions were fitted to a combination of a Gaussian

Solvation Dynamics of C153 in Supercritical Fluoroform

J. Phys. Chem. B, Vol. 110, No. 20, 2006 10125

Figure 5. Solute-solvent pair correlation functions at three different densities for six C153 atoms: the first peak is shown in detail. The groundstate results are shown in full lines, while the excited-state results are shown in dashed lines.

Figure 6. Five-site fluoroform: solute-solvent orientational correlations in the ground and in the excited state at four different densities. The internuclear solute vector is along the C3-C8 direction.

and two exponential functions:

1 F(t) ) a exp - ω2t2 + b1 exp(-t/τ1) + b2 exp(-t/τ2) 2 (15)

(

)

Equation 15 does not have a precise theoretical justification, but it is often employed to describe solvation response obtained from simulation because it has the correct long-time behavior and close to the correct behavior at short times.11 The best results for the fit were obtained by using the constraint a + b1 + b2 ) 1, and they are shown in the bottom panels of Figures 3 and 4.

Figure 7. Same as in Figure 6, but for the N13-C2 direction.

TABLE 5: Fits to Eq 15, Ground State a ω2 (ps-2) b1 τ1 (ps) b2 τ2 (ps)

0.13 Fc

0.3 Fc

0.59 Fc

1.5 Fc

2.0 Fc

0.52 28.4 0.30 0.47 0.18 3.52

0.68 23.6 0.34 0.16 0.18 2.04

0.57 26.7 0.28 0.32 0.15 1.56

0.54 30.3 0.36 0.40 0.10 1.20

0.51 33.2 0.34 0.28 0.15 0.91

The ground-state parameters are shown in Table 5. The decay is dominated by the Gaussian portion, even though the contribution to the total response gets smaller as the density increases. The contribution (coefficient b2 in eq 15) associated with the slower exponential (relaxation time τ2) is more or less constant, and it is smaller than the contribution (b1) due to the faster

10126 J. Phys. Chem. B, Vol. 110, No. 20, 2006

Ingrosso and Ladanyi

Figure 8. Five-site fluoroform: comparison between the equilibrium and the nonequilibrium solvation response at two different densities.

TABLE 6: Fits to Eq 15, Excited State a ω2 (ps-2) b1 τ1 (ps) b2 τ2 (ps)

0.13 Fc

0.3 Fc

0.59 Fc

1.5 Fc

2.0 Fc

0.58 38.7 0.20 1.39 0.22 30.6

0.45 39.2 0.20 0.44 0.20 15.0

0.45 28.1 0.35 0.28 0.20 5.38

0.52 27.6 0.28 0.20 0.20 2.45

0.48 34.2 0.33 0.25 0.19 1.14

exponential. τ2 becomes smaller as density increases, while the density dependence of τ1 is less regular. The excited-state parameters are reported in Table 6. Compared to the ground-state results, the relative contribution of each portion to the total response has a weaker density dependence. The relaxation time related to the slower exponential shows a strong density dependence, becoming much smaller at high density. The slower exponential dominates the diffusive portion of the decay, which results in a stronger density dependence of the excited-state results compared to the groundstate ones. In the ground-state response, the slower relaxation time does not show a dramatic change with density compared to the excited-state case, where it varies by 1 order of magnitude from the lowest to the highest density. The faster relaxation time decreases at intermediate densities and then increases slightly with increasing density. We present in Figure 8 a comparison between the equilibrium and the nonequilibrium response at two densities above and below the critical density. In both cases, the nonequilibrium response lies between the ground and the excited-state equilibrium ones, and the relaxation rates of all functions become more similar at higher densities. In the lower density case, the nonequilibrium response decays with a rate that is more similar to the equilibrium excited-state response in the longer time scale. This finding is in agreement with previous simulations,8,11,12 The comparison with experimental data of Kometani, Arzhantsev, and Maroncelli (KAM)15 for DCS in SC fluorform can be performed for F g 1.4Fc. We use the parameters in Tables 5 and 6 for equilibrium ground- and excited-state responses to carry out the comparison. We found that the largest contribution to the solvation response is Gaussian (about 50% of the overall decay for both equilibrium and nonequilibrium responses). On the other hand, experimental data were fit to a sum of two exponentials and the Gaussian function was not included. The Gaussian portion of the decay is needed in the fit of computed solvation response for a good description of the ultrafast time scale of solvation dynamics. Depending on the pulse width used in experiments,

Figure 9. The effect of the inclusion of induced-dipole interaction in the calculation of the solvation response of C153 in SC fluoroform at 1.5 Fc compared to experimental measurements15 of DCS in the same solvent.

the very short portion of the time evolution might be not detectable, generating discrepancies with simulations.26 In the experiments, the larger contribution to the overall response is due to the faster exponential (about 80%), which accounts for about 30-35% of the total response in the simulation. The value of the relaxation time is smaller in simulation than in experiments. Presumably the faster exponential in experiments is the sum of the Gaussian and the faster exponential in our fit. The relative contribution of the slower exponential to the overall decay in ref 15 is around 15%. In our case, the slower exponential accounts for 20% of the total response in the excited state and for about 10-15% of the total response in the ground state, but the relaxation times are smaller than the experimental ones. We note here that the solute employed in the measurements (DCS) was not the same as the one employed in our simulation, and we cannot exclude a solute dependence of the solvation response. Direct measurements of solvation times for C153 in SC fluoroform have been obtained at two densities below critical (0.4 and 0.8 Fc) at 323 K.34 In that paper, it was estimated that 10-20% of the relaxation occurs with a 60-70 ps relaxation time. However, the resolution of the experiment (30 ps response time) is too low to use those results for a quantitative comparison with simulation. When we introduced the contribution of the induced-dipole interaction to ∆E according to eq 6, for a value of the bulk density corresponding to 1.5 Fc, we found that the calculated result is in better agreement with experiment for DCS in SC fluoroform. The comparison is shown in Figure 9. In the time scale corresponding to diffusive motions, the calculated ground- and excited-state responses approach the long time scale experimental result better than the corresponding results obtained without introducing the induced-dipole interaction. Improvement is also observed in the time scale connected with nondiffusive solvent reorientation. It is worth noting that, because we employed a rigid potential for the solvent in order to evaluate polarizability effects, the feature related to intramolecular bending is not present in the solvation response. From the analysis of the solvation response, it appears that the five-site model is better capable of a realistic description of the solvent reorientational dynamics. Further insights into model dependence of fluoroform dynamics are obtained by the calculation of orientational TCFs. 4.2. Orientational Relaxation. We start by presenting a comparison between the two-site and the five-site models based

Solvation Dynamics of C153 in Supercritical Fluoroform

Figure 10. Reorientational dynamics in SC fluoroform: comparison between results obtained with the two-site model16 (dotted-and-dashed line) and the five-site model17 (full line) for the first rank TCFs (eq 7, l ) 1, and eq 9, respectively).

on calculations of the first rank orientational time correlation functions. To make this comparison at the same temperature for both models, we had to run the two-site simulations at a temperature (T ) 318 K) slightly higher than the one corresponding to 1.03 Tc (T ) 310 K in the two-site case). In the five-site model, in fact, this temperature is too close to Tc (see Table 1). As displayed in Figure 10 for two of the densities considered, both the single-molecule and the collective orientational time correlations, defined in eq 7, decay more rapidly in the twosite case in the librational portion of their time evolution and have a smaller amplitude when the diffusive decay sets in. No differences are to be expected in the time scale associated with purely inertial motions, corresponding to the ultrafast, initial decay displaying a Gaussian shape. In fact, as already mentioned, the two-site model was developed to reproduce the real moment of inertia for rotation of the fluoroform symmetry axis, and therefore, the initial curvature of the Gaussian function must be the same at the same temperature.35 The result of faster dynamics for the two-site model in the time scale associated with reorientational motions confirms the result obtained for the solvation response. In particular, in Figure 2, the dip below zero occurs on a time scale characteristic of the nondiffusive reorientational motion of a fluoroform molecule. The result of the interaction with the surroundings is a torque acting on the molecule, which influences the reorientation and results in what is generally referred to as libration.36,37 The too-fast reorientational dynamics obtained from the two-site model suggests that the torque on each librating fluoroform molecule is weaker than it should be. The difference in intermolecular torques does not account for all the differences between the solvation responses of the two models because, as already discussed, the solvation response function also depends on translational motions. However, because improved agreement with experiment is obtained with the five-site model, we conclude that a realistic description of the molecular shape is necessary to study solvation dynamics in near-critical fluoroform. From now on, the discussion will be focused on the reorientational dynamics results obtained with the five-site model. In Table 7, we report the correlation times for the relaxation of the single-molecule (first and second rank, eq 8) and of the collective dipole (first rank, eq 10) TCFs.

J. Phys. Chem. B, Vol. 110, No. 20, 2006 10127

Figure 11. Five-site fluoroform: second-rank orientational time correlation functions from eq 7, with l ) 2.

TABLE 7: Rotational Relaxation Times Fr 0.13 0.3 0.59 1.5 2.0

τ1 (ps)

τ1c (ps)

τ2 (ps)

1.25 0.76 0.65 0.68 0.87

1.42 1.04 0.84 0.83 1.22

0.61 0.44 0.28 0.28 0.33

Both correlation times are characterized by a higher value at low densities by very small changes in the intermediate region and by an increased value for the liquidlike density. The single-molecule reorientational dynamics is very similar, even though slightly faster, than collective dynamics. The latter is more influenced by changing from intermediate to liquidlike densities because the change in the correlation time is larger than the single-molecule case. We turn now to consider rotational dynamics, as described by Ψ2(t) defined in eq 7, with l ) 2. We shall use the results of ref 16 for a comparison with the two-site model. The calculated Ψ2(t) functions that we obtained are reported in Figure 11. In ref 22, rotational correlation times are obtained from polarized and depolarized Raman scattering measurements. What is observed from experiments is that the shape of the measured functions shows strong density dependence, as in the case of the calculated ones displayed in Figure 11. For densities larger than the critical density, the experimentally measured decay behaves as expected for a liquidlike system with high collision frequency.38,39 The experimental behavior is reproduced by the calculation. At low densities, the TCF resembles that of a free-rotor, except at long times, at which it decays to zero instead of reaching a finite limiting value. We report in Figure 12 the rotational correlation times obtained for the rank two functions according to eq 8. The lowdensity limit of τ2 is infinite, given that the molecule becomes a free rotor.38 We see in fact a large increase in the calculated time for the smallest density. As for the liquidlike density, the correlation time increases due to the presence of increasing friction. This behavior is similar to that observed for the first rank, single-molecule and for the collective dipole TCFs. The agreement with the experimental data is satisfactory because the shape of the density dependence is reproduced, and quantitatively reasonable, even though our results at densities higher than 0.3 Fc give τ2 values that are 16% smaller than the experimental ones. Underestimated rotational times were reported in ref 16 for the two-site model (underestimated by 25%), and we show here that a five-site model introduces an improvement for describing this property as well. The experimental data

10128 J. Phys. Chem. B, Vol. 110, No. 20, 2006

Figure 12. A comparison of the second-rank single-molecule rotational correlation times obtained from MD simulations of five-site fluoroform (squares) and from experimental Raman scattering measurements22 (circles).

were calculated through direct integration (as defined in eq 8) by assuming exponential decay of the TCFs at t larger than 1.5 ps, which could introduce some discrepancy with our treatment (our TCFs were calculated up to 40 ps). Nonetheless, the rotational times provided in ref 22 were shown to be in good agreement with other experimental work at higher density (larger than 1.5 Fc) based on NMR measurements.40 Smaller rotational times for the two-site model were explained in terms of an underestimated rotational volume due to the use of a very simple geometrical model for fluoroform.16 However, it is worth observing that there are effects related to rotational friction, which both models neglect, due to their nonpolarizable nature. Polarization effects were shown to be necessary to have a correct picture for intermolecular interactions in fluoroform.41 In ref 5, improved agreement between simulation and experiment was obtained when induction was taken into account to calculate the solvatochromic shifts from which local density augmentation was evaluated. Because the two-site model has an enhanced dipole16 (the molecular dipole moment was equal to 2.2 D, while the experimental value is 1.649 D42), it includes induction effects at a mean-field level. The five-site model does not have any mean-field correction, and the molecular dipole moment is equal to 1.67 D. 5. Concluding Remarks In this paper, we examined the dynamics of solvation in supercritical fluoroform. The investigation was performed by studying the solvation response function, describing the orientational rearrangement of the solvent molecules after a perturbation in the solute electronic structure and the reorientational dynamics of the pure solvent. Molecular dynamics simulations were employed to calculate the properties of the systems of interest. The solute, C153, was described through an all-atom potential model, while we compared results from two different models for the solvent and for the pure SC fluid: a two-site and an atomistic model. The differences obtained in the two cases are large, especially for the nondiffusive portion of the solvent response, which accounts for librational motion of solvent molecules. The discrepancy between the results obtained with different description of the solvent molecule suggests that the molecular shape plays an important role in obtaining an accurate picture of dynamics. The expectation, based on earlier computer simulations of solvation dynamics in the SC regime,8,11,12 that changes in the effective local density would lead to deviations from the linear

Ingrosso and Ladanyi response regime, was confirmed for the five-site model. However, because of smaller changes in the solvation shell structure in the present case, the deviations were smaller than those found in the other model SC solute-solvent systems. Our results indicate that the five-site model provides an improved description of reorientational dynamics compared to the two-site model. This improved description impacts solvation dynamics, as the comparison with experiments for DCS in SC fluoroform shows. We adopted a methodology based on a perturbative approach to include the interaction between the solute charge distribution and the induced dipoles on the solvent in the calculation of the time-dependent solvation energy. The result thus obtained was in better agreement with experiment, suggesting that induction effects need to be taken into account for a good description of interactions in the SC regime. To provide an explanation for the disagreement between the two-site and the five-site model results for solvation dynamics, we studied reorientational dynamics as described by the singlemolecule and collective dipole orientational TCFs. By comparing the results for the first rank functions at two different densities, what we learn is that the orientational relaxation dynamics is too fast in the two-site model compared to the five-site model. In the paper presenting the model,16 the authors provided an explanation for this by relating the faster rotational relaxation to a smaller rotational volume compared to the molecular one. This, in addition to a more realistic molecular shape, might be the reason for a better description provided by the five-site model. A model with unequal Lennard-Jones parameters for the two sites might provide a more realistic two-site representation of SC fluoroform. A further test for the five-site model was provided by comparing our calculated rotational times with experiments,22 which showed fairly good agreement. We can therefore conclude that accurate calculations of timedependent properties in the SC regime for fluoroform require a realistic representation of the molecular shape. Additional improvement can be obtained by including nonadditive effects in the intermolecular interaction by means of polarizable potentials. Acknowledgment. This work was supported by the NSF through grant CHE9981539. We are grateful to Prof. Mark Maroncelli and to Dr. Noritsugu Kometani for providing us their experimental results on solvation dynamics in supercritical fluoroform. We also thank Prof. Maroncelli for carefully reading the manuscript and providing us with his comments and suggestions. References and Notes (1) Tucker, S. C. Chem. ReV. 1999, 99, 391. (2) Kajimoto, O. Chem. ReV. 1999, 99, 355. (3) Sunol, A. K.; Sunol, S. S. In Handbook of SolVents. Substitution of SolVents by Safer Products and Processes; Wypych, G., Ed.; William Andrew and Chem Tec Publishing: Toronto, Canada, New York, 2001. (4) Biswas, R.; Lewis, J. E.; Maroncelli, M. Chem. Phys. Lett. 1999, 310, 485. (5) Ingrosso, F.; Ladanyi, B. M.; Mennucci, B.; Scalmani, G. J. Phys. Chem. B 2006, 110, 4953. The vertical axis labels of the top and bottom panels in Figures 6 and 7 of this paper were inadvertently interchanged. The top panel should be Feff/F, and the bottom panel label ∆Feff,r. (6) Kapko, V.; Egorov, S. A. Chem. Phys. Lett. 2005, 402, 258. (7) Kapko, V.; Egorov, S. A. J. Chem. Phys. 2004, 121, 11145. (8) Egorov, S. A. J. Chem. Phys. 2004, 121, 6948. (9) Egorov, S. A. J. Chem. Phys. 2000, 113, 1950. (10) Egorov, S. A. J. Chem. Phys. 2002, 116, 2004. (11) Ladanyi, B. M.; Nugent, S. J. Chem. Phys. 2006, 124, 044505. (12) Re, M.; Laria, D. J. Phys. Chem. B 1997, 101, 10494. (13) Martins, M. M.; Stassen, H. J. Chem. Phys. 2003, 118, 5558. (14) Graf, P.; Nitzan, A. Chem. Phys. 1998, 235, 297.

Solvation Dynamics of C153 in Supercritical Fluoroform (15) Kometani, N.; Arzhantsev, S.; Maroncelli, M. J. Phys. Chem. A 2006, 110, 3405. (16) Song, W.; Patel, N.; Maroncelli, M. J. Phys. Chem. B 2002, 106, 8783. (17) Potter, S. C.; Tildesley, D. J.; Burgess, A. N.; Rogers, S. C. Mol. Phys. 1997, 92, 825. (18) Horng, M. L.; Gardecki, J. A.; Papazyan, A.; Maroncelli, M. J. Phys. Chem. 1995, 99, 17311. (19) Kumar, P. V.; Maroncelli, M. J. Chem. Phys. 1995, 103, 3038. (20) Jimenez, R.; Fleming, G. R.; Kumar, P. V.; Maroncelli, M. Nature 1994, 369, 471. (21) Fleming, G. R.; Cho, M. Annu. ReV. Phys. Chem. 1996, 47, 109. (22) Okazaki, S.; Matsumoto, M.; Okada, I.; Maeda, K.; Kataoka, Y. J. Chem. Phys. 1995, 103, 8594. (23) Maroncelli, M. J. Chem. Phys. 1991, 94, 2084. (24) Carter, E. A.; Hynes, J. T. J. Chem. Phys. 1991, 94, 5961. (25) Nugent, S.; Ladanyi, B. M. J. Chem. Phys. 2004, 120, 874. (26) Ingrosso, F.; Ladanyi, B. M.; Mennucci, B.; Elola, M. D.; Tomasi, J. J. Phys. Chem. B 2005, 109, 3553. (27) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: London, 1976. (28) Smith, W.; Forester, T. R. DL_POLY_2; CCP5; Daresbury: Warrington, England, 1999.

J. Phys. Chem. B, Vol. 110, No. 20, 2006 10129 (29) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. (30) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (31) Ladanyi, B. M.; Maroncelli, M. J. Chem. Phys. 1998, 109, 3204. (32) Ladanyi, B. M. In Theoretical Methods in Condensed Phase Chemistry; Schwartz, S. D., Ed.; Kluwer: Dordrecht, The Netherlands, 2000; Chapter 7. (33) Stell, G.; Patey, G. N.; Høye, J. S. AdV. Chem. Phys. 1981, 48, 183. (34) Kimura, Y.; Hirota, N. J. Chem. Phys. 1999, 111, 5474. (35) St. Pierre, A. G.; Steele, W. A. Mol. Phys. 1981, 43, 123. (36) Lynden-Bell, R. M.; Steele, W. A. J. Phys. Chem. 1984, 88, 6514. (37) Steele, W. A. J. Chem. Phys. 1963, 38, 2404. (38) St. Pierre, A. G.; Steele, W. A. J. Chem. Phys. 1972, 57, 4638. (39) Gordon, R. G. J. Chem. Phys. 1966, 44, 1830. (40) Lang, E. W.; Prielmeier, F. X.; Radkowitsch, H.; Ludemann, H. D. Ber. Bunsen.-Ges. Phys. Chem. 1987, 91, 1025. (41) Hloucha, M.; Deiters, U. K. Fluid Phase Equilib. 1998, 149, 41. (42) Gershel, A.; Dimicoli, I.; Haffre, J.; Riou, A. Mol. Phys. 1976, 32, 679. (43) Penoncello, S.; Shan, Z.; Jacobsen, R. ASHRAE Trans. 2000, 106, 739.