Smoluchowski Equation Description of Solute Diffusion Dynamics and

18060

J. Phys. Chem. C 2007, 111, 18060-18072

Smoluchowski Equation Description of Solute Diffusion Dynamics and Time-Dependent Fluorescence in Nanoconfined Solvents Xiaobing Feng and Ward H. Thompson* Department of Chemistry, UniVersity of Kansas, Lawrence, Kansas 66045 ReceiVed: June 11, 2007; In Final Form: August 23, 2007

The diffusion dynamics and time-dependent fluorescence of a model dye molecule dissolved in CH3I, CH3CN, and CH3OH solvents confined in nanoscale cavities and pores of 20-30 Å diameter have been modeled using a Smoluchowski equation approach. The effects of the confining framework geometry (spherical or cylindrical) and size have been investigated. The two-degree-of-freedom free-energy surfaces for the dye molecule ground and excited states are constructed from one-dimensional free-energy curves in the solute position and model force constants for a collective solvent coordinate. Multiple time scales are found in the solute diffusion after excitation and in the time-dependent fluorescence. The origins of these time scales can be understood in terms of solvent reorganization and dye molecule diffusion. The results indicate that the Smoluchowski equation approach has promise for describing time-dependent fluorescence dynamics, and ultimately reaction dynamics, in a wide variety of confined solvent systems.

1. Introduction There is currently significant interest in nanostructured materials. In particular, microporous and mesoporous materials, which can confine liquids in nanoscale pores, have been the focus of both synthetic and characterization efforts. The former indicate the promise of designing such materials for a variety of applications,1-3 and the latter are demonstrating the significant effects on chemistry that are associated with nanoscale confinement of a solvent. In particular, there have been many studies of time-dependent fluorescence (TDF).4-12 In these experiments, a dye molecule with an electronic charge-transfer transition is excited and the resulting fluorescence is monitored as a function of time. The changes in fluorescence energy with time provide information about the reorganization of the solvent around the dye molecule in response to a change in its dipole moment upon excitation. The solvation dynamics probed are closely related to those involved in other charge-transfer processes such as electron- or proton-transfer reactions. Moreover, upon nanoconfinement the solvation dynamics are often changed dramatically: in bulk solvents the solvent reorganization is typically complete after a few picoseconds, whereas in nanoconfined solvents it involves multiple time scales and can have slow components on the order of nanoseconds. Understanding how confinement and the properties of the confining framework (e.g., size, dimensionality, flexibility, surface functionality) affect the solvation dynamics is important in designing materials for specific purposes including as reaction media. A number of models have been proposed to describe the multiexponential, long-time decay observed in time-dependent fluorescence studies of nanoconfined solvents. One of the most prevalent models for describing the behavior of confined solvents is a two-state model in which a property is assumed to originate from the superposition of the properties of solvent molecules in different regions.4,13 Specifically, the picture is that the solvent molecules in the “surface” layer at the solventconfining framework interface respond on a slower time scale * Corresponding author. E-mail: [email protected].

than the “interior” solvent molecules away from the interface. This model is sometimes used to interpret TDF measurements, though some questions about its applicability to other observables have been raised recently,14,15 Another model, based on exchange of surface and interior water molecules, has been proposed by Bagchi and co-workers.16 Finally, one of us recently proposed that long time scales may arise in TDF experiments in nanoconfined solvents because of diffusion of the dye molecule after excitation.17 This is a result of the positiondependent solvent polarity: the nanoconfined solvent is effectively less polar near the interface with the confining framework and increases in polarity with the distance from that surface.18 This leads to electronic state-dependent position distributions for the dye molecule, which typically has a small dipole moment in the ground state, and hence sits near the interface, and a large dipole moment in the excited state, and so is found away from the interface.18,19 It is not immediately clear which of these models, if any, applies for a given nanoconfined solvent. It may be that the multiexponential decays of the TDF Stokes shift have physical origins that depend on the nature of the confining framework, solvent, and even possibly the dye molecule. Although molecular dynamics (MD) simulation provides a particularly appealing avenue for exploring the connection between the system properties and the characteristics of the TDF signal, full nonequilibrium simulations can be quite time-consuming. In nanoconfined solvents, the effort is compounded by the heterogeneity of the system and the long-time dynamics, which together mean that many lengthy nonequilibrium trajectories are required. This provides an impetus for developing other ways of screening the TDF dynamics rapidly for a large number of nanoconfined solvent systems. To that end, we describe in this paper the application of the Smoluchowski equation approach to the description of TDF in solvents confined in spherical cavities and cylindrical pores. In addition to providing a rapid approach to simulating the time-dependent Stokes shift and solute diffusive motion, the Smoluchowski equation has a number of other advantages.

10.1021/jp074516h CCC: $37.00 © 2007 American Chemical Society Published on Web 11/13/2007

Smoluchowski Equation Description

J. Phys. Chem. C, Vol. 111, No. 49, 2007 18061

Specifically, the entire time-dependent fluorescence spectrum is obtained as a function of time; the sampling in nonequilibrium MD simulations is rarely sufficient to reconstruct the entire spectrum. The solute position distribution function is also calculated as a function of time in the approach, allowing insight into the mechanism(s) of diffusion (vide infra). Furthermore, it is relatively easy to change the geometry of the confining framework within the Smoluchowski equation so that spherical, cylindrical, and slit pore geometries can be compared straightforwardly. However, the approach does require knowledge of the potential of mean force for the solute, in general as a function of both the solute position and a collective solvent coordinate. In this work, we use a combination of calculated freeenergy curves and reasonably constructed model free-energy surfaces. The organization of the remainder of this paper is as follows. The Smoluchowski equation formalism for spherical and cylindrical confining frameworks is presented in Section 2.1. The free-energy surfaces, diffusion coefficients, and other parameters defining the systems studied in this work are presented in Section 2.2 along with the details of the simulations. The results for the time-dependent fluorescence of a model dye molecule dissolved in CH3I, CH3CN, and CH3OH solvents confined in spherical cavities and cylindrical pores of 10 and 15 Å radius are presented and discussed in Section 3. The effect of cavity/pore size, free-energy surface characteristics, diffusion coefficients, and geometry (spherical vs cylindrical) are all addressed. Finally, some concluding remarks are offered in Section 4. 2. Theory and Models 2.1. Smoluchowski Equation. The Smoluchowski equation is the high-friction limit of the Fokker-Planck equation.20 Specifically, it is obtained by assuming that the momentum relaxation is instantaneous, and it is valid on diffusive time scales.21 The Smoluchowski equation for a solute molecule in a nanoconfined solvent can be written as22

∂n(r, t) ) ∇[D(r)(∇n(r, t) + βn(r, t)∇A(r))] ∂t

(1)

where n(r, t) is the probability density at position r and time t, D(r) is the position-dependent diffusion coefficient, β ) 1/kBT, T is the temperature, and A(r) is the free energy of the system with the center of mass of the solute at position r. In this paper, we consider spherical cavities and cylindrical pores. In these cases, the free energy depends only on the radial position, so A(r) can be reduced to a one-dimensional free-energy curve, A(r), where r is the distance of the solute from the cavity or pore center. When a solute chromophore in a bulk solvent is excited electronically the solvent reorganizes to achieve favorable solvation of the new charge distribution. This process, which is solvent-dependent, typically occurs on the time scale of a few picoseconds or less. To include this effect, which is not present in eq 1, the solvent coordinate, ∆E, is introduced as an additional argument of the probability density. The time evolution of the probability density as a function of position r and the solvent coordinate would enable us to study solute radial diffusion and fluoresence at the same time. The solvation (excluding the initial ultrafast inertial component) is mainly from the diffusive rotational motion of the solvent molecules; thus, the solvation dynamics can be described by a diffusive motion, as many workers have done previously.23-28 The corresponding

equation for the joint probability density is obtained by generalizing eq 1:

∂n(r, ∆E, t) ) ∇[D(r)(∇n(r, ∆E, t) + ∂t βn(r, ∆E, t)∇U(r, ∆E))] + ∂U(r, ∆E) ∂ ∂n(r, ∆E, t) (2) DE + βn(r, ∆E, t) ∂∆E ∂∆E ∂∆E

(

)

Here, DE is a diffusion constant dictating the time scale on which the solvent reorganization happens, U(r, ∆E) is the multidimensional free-energy surface, which can be reduced to a function of two coordinates, U(r, ∆E), in the case of a spherical cavity or cylindrical pore. The collective solvent coordinate, ∆E, can be defined in many ways; in this work, ∆E ) Ve(Q) - Vg(Q), where Vg(e)(Q) is the energy of the system with the solute in ground (excited) state, and Q represents the coordinates of all atoms in the system save the solute center-of-mass position. In modeling a TDF experiment, the solute is assumed to initially follow a Boltzmann distribution determined by the ground-state free-energy surface; the equilibrium (long-time) distribution of the solute after electronic excitation should also be a Boltzmann distribution determined by the excited-state freeenergy surface. It is easy to verify that eq 2 with U(r, ∆E) the excited-state free-energy surface has the correct equilibrium distribution, that is, a Boltzmann distribution as t f ∞. For spherical cavities and cylindrical pores, the free energies and initial distributions have spherical or cylindrical symmetry that is preserved during the time evolution. Thus, as noted above, the Smoluchowski equations can be reduced to radial versions. Therefore, a new distribution, f(r, ∆E, t), is introduced as f(r, ∆E, t) ) 4πr 2n(r, ∆E, t) for a spherical cavity and f(r, ∆E, t) ) 2πrn(r, ∆E, t) for a cylindrical pore. The corresponding Smoluchowski equation can be written for both cases as

{

} }

∂f(r, ∆E, t) ∂ ∂ ) D(r) e-βV(r,∆E) [eβV(r,∆E) f(r, ∆E, t)] + ∂t ∂r ∂r ∂ ∂ [eβV(r,∆E) f(r, ∆E, t)] (3) e-βV(r,∆E) DE ∂∆E ∂∆E

{

where V(r, ∆E) ) U(r, ∆E) - 2kBT ln r for a spherical cavity, and V(r, ∆E) ) U(r, ∆E) - kBT ln r for a cylindrical pore.29 Note that the normalization condition for both distributions is ∫∫f(r, ∆E, t)drd∆E ) 1. It is often observed that the free-energy curve in ∆E for fixed r is close to parabolic;18,30-32 therefore, it is a good approximation to write Uj(r, ∆E), the free-energy surface for state j, as

1 Uj(r, ∆E) ) Uj0(r) + kj(r)[∆E - Ej0(r)]2 2

(4)

where kj(r) is the position-dependent force constant, Ej0(r) is the position-dependent minimum in the solvent coordinate, and Uj0(r) is the minimum free energy for fixed r. It may appear that specifying the two-dimensional free-energy surfaces for both ground and excited states requires the knowledge of six functions, that is, Uj0(r), kj(r), and Ej0(r) for j ) g (ground state) and j ) e (excited state). However, as is shown below, with calculated one-dimensional free-energy curves, Aj(r), as inputs, Uj0(r) and Ej0(r) are determined by the knowledge of kj(r) and the requirement33 that ke(r) ) kg(r). The free energy, A(r), is defined as34

18062 J. Phys. Chem. C, Vol. 111, No. 49, 2007

1 A(r) ) - ln β

∫ e-βV(r,Q)dQ ∫ e-βV(r′,Q′)dr′dQ′

Feng and Thompson

(5)

where V is the system potential energy, r and r′ denote solute positions, and Q (Q′) represents the remaining coordinates of the system. Similarly, U(r, ∆E) is given by

1 U(r, ∆E) ) - ln β

∫ e-βV(r,Q)δ[S(Q) - ∆E]dQ ∫ e-βV(r′,Q′)dr′dQ′

(6)

where S(Q) is the difference in solute-solvent interaction energies with the solute in the excited state and with the solute in the ground state, that is, the instantaneous value of the solvent coordinate. Using eqs 5 and 6, one can obtain the following relation, ∫e-βU(r,∆E)d∆E ) e-βA(r). On substituting eq 4 into the above relation, one then has

Uj0(r) ) Aj(r) +

1 2π ln 2β βkj(r)

(7)

Moreover, Tachiya has shown that in a bulk solvent the free energy of the excited state differs from the free energy of the ground state by ∆E.31,33 An analogous relationship applies to confined solvents; it is (see Appendix)

Ue(r, ∆E) ) ∆E + Ug(r, ∆E) + Ug - Ue

(8)

Here, U h e and U h g are the total free energies of the system averaged over the solute position and solvent coordinate in the excited and ground electronic states, respectively. Using eqs 4, 7, and 8, the following relationships can be obtained

Ue(r, ∆E) ) Ae(r) +

2π 1 1 ln + k(r)[∆E - Ee0(r)]2 (9) 2β βk(r) 2

Ug(r, ∆E) ) Ag(r) +

1 2π + ln 2β βk(r) 1 k(r)[∆E - Eg0(r)]2 (10) 2

Ee0(r) ) Eg0(r) -

1 k(r)

Eg0(r) ) Ae(r) - Ag(r) +

1 2k(r)

(11) (12)

It can be seen that the force constants are the same for the excited and ground states, kg(r) ) ke(r) ≡ k(r), which is a direct consequence of eq 8. The locations of the minima for the ground-state free-energy surface are determined by

∂Ag(r) 1 ∂k(r) )0 ∂r 2βk(r) ∂r

(13)

An analogous equation applies to the excited state. If the force constant is independent of radial position, then the minima are exactly the same as those in the one-dimensional free-energy curves. For most reasonable choices of k(r), the changes in the radial positions of the minima are small. Therefore, the global minima of the free-energy surfaces have nearly the same radial coordinates as for the corresponding one-dimensional freeenergy curves.

Figure 1. Position-dependent force constants generated from a quadratic function, k(r) ) a + br + cr 2. a ) 0.119 (kcal/mol)-1 is used for all cases. The other parameters are: b ) -0.0054, c ) 0.0011, for case I (black); b ) 0, c ) 0, for case II (red); b ) -0.0054, c ) 0.0031, for case III (blue); b ) -0.016, c ) 0.0011, for case IV (green). b and c are in units of (kcal/mol)-1 Å-1 and (kcal/mol)-1 Å-2, respectively. The k(r) shown here are for R ) 10 Å.

2.2. Simulation Details. 2.2.1. Free-Energy Surfaces. Three solvents, CH3I, CH3CN, and CH3OH, are considered here, and the Ag(r) and Ae(r) for the solute confined in these solvents calculated previously by Mitchell-Koch and Thompson35 using thermodynamic integration are used in this work. The solute is modeled as a diatomic molecule, with a bond length equal to 3 Å. The two atoms of the model solute differ only in charge. In its ground state, the two atoms have charges of (0.1 |e|, which gives a small dipole moment of 1.44 D. When the solute is electronically excited, the dipole moment increases to 7.1 D because of atomic charges of (0.5 |e| charges. The LennardJones parameters are ) 1.6637 kcal/mol and σ ) 3.5 Å. The mass of each solute atom is 40 amu. The CH3 groups in the solvent molecules are treated as united atoms. Intermolecular potentials consist of Lennard-Jones and Coulombic interactions. The nanoconfined solvent density is taken to be 90% of the corresponding bulk solvent density.35 Solvents and the solute interact with the cavity wall only through an integrated LennardJones potential;36,37 therefore, the cavity is hydrophobic. The Lennard-Jones parameters for the cavity wall are wall ) 3 kcal/ mol and σwall ) 2.5 Å. Further details about the models for the solute, the solvents, and the cavity can be found in refs 17, 18, and 35. As noted above, the two-dimensional free-energy surfaces, that is, Ue(r, ∆E) and Ug(r, ∆E), can be constructed from a model k(r) and the one-dimensional free-energy curves, Ae(r) and Ag(r). Limited results for k(r) have been obtained for some systems in previous Monte Carlo simulations.18 To investigate the effect of the position-dependence of the force constant, we consider four cases for k(r), shown in Figure 1. The k(r) of Case I is based on fitting previous Monte Carlo simulation data18 to a quadratic function; k(r) of Case II is a constant with a value close to the k(r) of Case I. The k(r) for Cases III and IV allows us to study the effects of large, qualitatively different, variations of the force constants that may represent other systems with stronger surface interactions. The two-dimensional free-energy surfaces, Ug(r, ∆E) and Ue(r, ∆E), constructed from one-dimensional free-energy curves, Ag(r) and Ae(r), and the four k(r) cases shown in Figure 1 are plotted in Figure 2. As is expected, the global minima of the free energies of the ground state lie near the cavity wall, whereas those of excited-state lie near the cavity center. The free-energy surfaces diverge both near the cavity center and wall because



Figure 2. Contour plots of two-dimensional free-energy surfaces of ground (blue) and excited (red) states. The solvent is CH3I confined in a spherical cavity with radius equal to 10 Å. The free-energy surfaces are generated by using the four different k(r) in Figure 1. The free-energy surfaces corresponding to Cases I-IV of k(r) are shown in panels a-d, respectively. The x axis is the solvent coordinate in kcal/mol. The y axis is the distance from cavity center in Å. The spacing is 0.3 kcal/mol.

Figure 3. Contour plots of initial (blue) and equilibrium (red) distributions corresponding to the free-energy surfaces shown in Figure 2. The x axis is the solvent coordinate in kcal/mol. The y axis is the distance from cavity center in Å. The solvent is CH3I; the cavity is spherical with a radius equal to 10 Å.

the free-energy curves are diverging at these places. For nanoscale cavities, the solvent molecules are not distributed uniformly but rather are packed in layers. The free energy for the solute near the center of the cavity is thus large because the solvent layering is incommensurate with having a molecule

there. The corresponding initial and equilibrium distributions are shown in Figure 3. 2.2.2. Solution of the Smoluchowski Equation. The Smoluchowski equations for f1D(r, t), which is the radial distribution function corresponding to the radial version of eq 1, and


Feng and Thompson

Figure 4. Comparison of the results of the one-dimensional Smoluchowski equation (eq 1) with that of the two-dimensional equation (eq 2) with different DE. The radial distribution, f(r, t), at t ) 100 ps is plotted in the left panel, and the average radial displacement, 〈∆d(t)〉, (see the text) is plotted in the right panel. Results are shown for DE ) 2.735 × 10-5 (kcal/mol)fs-1 (black), 2.735 × 10-4 (kcal/mol)fs-1 (red), 2.735 × 10-3 (kcal/mol)fs-1 (blue), and the one-dimensional equation (green), respectively. It can be seen that the one-dimensional solution is the large DE limit of solutions of the two-dimensional equation. The solvent is CH3I; the cavity is spherical with radius equal to 10 Å.

f(r, ∆E, t) are solved numerically using a finite difference method.38 Because the free-energy surfaces diverge at both large and small r, a natural boundary condition is used; that is, the distribution is set to zero at the boundaries specified by r ) 0 Å, r ) R - 2.5 Å, ∆E ) -20 kcal/mol, and ∆E ) 25 kcal/ mol, where R (10 or 15 Å) is the cavity or pore radius. The free energies at these boundaries are high enough to validate such an approach. Note also that before the solute is excited electronically its distribution is the Boltzmann distribution corresponding to the ground-state free-energy surface; this is the initial condition for the equations. The grid spacings are 0.024 Å in r and 0.15 kcal/mol in ∆E, and the time step is 0.1 fs; these choices give accurate results (doubling the grid points resulted in little change in the calculated quantities). All simulations are run to 1 ns. For the solvent CH3CN, which shows longer time scales than CH3I and CH3OH, the time scales extracted from a 2 ns simulation are very close to the ones from a 1 ns run. 2.2.3. SolVent Reorganization Time-Scale. The effective diffusion coefficient DE is related to the time scale of solvent reorganization, so it is important to consider its value for the systems of interest in this work and whether a separation of time scales between solvent reorganization and solute diffusion exists. A measure of this is obtained by comparing the radial distribution, f1D(r, t), calculated from the one-dimensional Smoluchowski equation, eq 1, to the fr(r, t) calculated from integrating over ∆E the distribution f(r, ∆E, t) resulting from the two-dimensional equation, eq 2

fr(r, t) )

∫-∞∞ f(r, ∆E, t)d∆E

(14)

In general, the radial distributions from eqs 1 and 2 are different, with the former being the large DE limit of the latter. A comparison of f(r, t) obtained from the one-dimensional and two-dimensional equations is shown in Figure 4a. The average change in solute distance from the cavity wall after excitation, 〈∆d(t)〉 ) ∫f(r, t ) 0)rdr - ∫f(r, t)rdr, calculated from the f(r, t) of Figure 4a are shown in Figure 4b (f(r, t) represents either f1d(r, t) or fr(r, t)). It can be seen from Figure 4 that when DE is increased the results from eq 2 approach the ones from eq 1. Better agreement with the one-dimensional results is obtained for DE ) 2.735 × 10-2 (kcal/mol) fs-1, the results of

which are not shown in Figure 4 because they overlap with the one-dimensional results. A longer solvent reorganization time (smaller DE) tends to slow down the diffusive motion of the solute because the solvent reorganization lowers the barriers the solute must diffuse across. In the rest of the work, DE is taken to be 2.735 × 10-3 (kcal/ mol) fs-1 to reproduce the intermediate time scale (∼1.7 ps) associated with solvent reorganization in previous nonequilibrium TDF simulations.17 Thus, the systems we consider have a significant separation of time scales for solvent reorganization and solute diffusion; the radial distribution and related quantities corresponding to this choice of DE are close to the ones from eq 1. Nonetheless, all of the results reported henceforth are calculated using eq 2 to avoid a needless approximation. 2.2.4. Position-Dependent Diffusion Coefficients. Another parameter appearing in eq 2 is the solute diffusion coefficient, D(r), which is, in principle, dependent on the solute position in the cavity or pore. To determine reasonable values for D(r), we used two approaches for the CH3I and CH3CN solvents in spherical cavities: (1) the results of the one-dimensional Smoluchowski equation were fit to nonequilibrium MD data, and (2) the diffusion coefficient was calculated from the Einstein relation37,39 using equilibrium molecular dynamics simulations. In the case of the latter, the position dependence, D(r), was estimated by calculating the root-mean-square (RMS) displacements separately for the solute within two regions of the cavity (defined by the radial distance from the cavity center). Using the Einstein relation,39 we calculated the diffusion coefficient in a bin from the slope obtained from a linear fitting to the RMS displacements (over the time range of 8-15 ps). Using the first approach for the solute dissolved in CH3I confined in hydrophobic spherical cavities, fitting the Smoluchowski equation results to the average radial displacement of the solute after excitation found in previous nonequilibrium MD simulations,17 gave D ) 4.5 × 10-6 cm2/s for R ) 10 Å and 9.0 × 10-6 cm2/s for R ) 15 Å. To determine if these values are reasonable and obtain an estimate for the positiondependence of D, we used the second approach. Equilibrium MD simulations did give results that were generally consistent, if not in exact agreement, with the fitting, 4.38 ( 0.10 × 10-6 cm2/s for R ) 10 Å and 7.48 ( 0.06 × 10-6 cm2/s for R ) 15 Å. The value for the smaller cavity results from an average of



TABLE 1: Parameters for Determining D(r) of CH3I and CH3CNa solvent

R (Å)

Case

D0 (cm2/s)

∆D (cm2/s)

CH3I CH3I CH3I CH3I CH3I CH3I CH3CN CH3CN CH3CN CH3CN CH3CN CH3CN

10 10 10 15 15 15 10 10 10 15 15 15

I-IV V VI I-IV V VI I-IV V VI I-IV V VI

4.5 × 10-6 4.5 × 10-6 4.5 × 10-6 9.0 × 10-6 9.0 × 10-6 9.0 × 10-6 2.64 × 10-6 2.64 × 10-6 2.64 × 10-6 5.8 × 10-6 5.8 × 10-6 5.8 × 10-6

0 1 × 10-6 5 × 10-7 0 1.5 × 10-6 7.5 × 10-7 0 8 × 10-7 4 × 10-7 0 1.2 × 10-6 6 × 10-7

R (1/a0)b 0.01 1 0.01 1 0.01 1 0.01 1

a The other parameters for Cases V and VI are the same as those for Case I. b a0 is the Bohr radius.

the diffusion coefficient of 3.57 ( 0.03 × 10-6 cm2/s for the solute near the cavity center (0 < r < 3.9Å) and the diffusion coefficient of 4.60 ( 0.07 × 10-6 cm2/s for the solute nearer the wall (r > 3.9 Å).40 For the larger cavity, diffusion coefficients of 6.56 ( 0.05 × 10-6 cm2/s for 0 < r < 7 Å, and 7.98 ( 0.06 × 10-6 cm2/s for 7 < r < 15 Å were obtained. An analogous analysis was used for CH3CN; however, no nonequilibrium MD simulations were available for R ) 10 Å. Thus, a diffusion coefficient of 3.1 × 10-6 cm2/s was estimated by extrapolating from data for larger cavities. The position dependence of the diffusion coefficient for R ) 10 Å was taken from equilibrium MD simulations of a 12 Å radius cavity, which gave a difference of 8 × 10-7 cm2/s between D near the cavity center and near the cavity wall. The diffusion coefficient of the solute in CH3CN confined in a 15 Å spherical cavity is directly calculated to be 6.8 ( 0.2 × 10-6 cm2/s, whereas a value of 5.8 × 10-6 cm2/s is found from fitting the nonequilibrium MD data. The calculated diffusion coefficient results from the average of 6.1 ( 0.3 × 10-6 cm2/s for r < 7 Å and 7.3 ( 0.3 × 10-6 cm2/s for 7 e r e 15 Å. The values from the fitting procedure were used with the relative magnitude of the position dependence from the equilibrium MD simulations to construct three models for D(r) of the form

D(r) ) D0 +

∆D tanh[R(r - R0)] tanh(RR/2)

(15)

Here R is the cavity radius, R0 ) (R - σwall)/2, D0 is the value of D(r) at r ) R0, ∆D measures the variation of D(r), and R controls the slope of D(r) at r ) R0. The model parameters are given for CH3I and CH3CN in Table 1; the diffusion coefficients obtained from fitting to nonequilibrium MD simulation results were used for D0 while the equilibrium MD simulations were used to estimate ∆D. Note that R ) 0.01 corresponds to essentially a linear change in D, whereas R ) 1 represents a smoothed step function. For CH3OH, no diffusion coefficients are available, so the values of CH3CN are used for this solvent. Therefore, the difference between time scales of the two solvents in this work comes completely from the difference in the free energies. Note that for each solvent and each geometry six cases are studied: In Cases I-IV, the corresponding force constants in Figure 1 are used along with a constant diffusion coefficient; in Cases V and VI, two position-dependent diffusion coefficients are used and the force constant of Case I in Figure 1 is used for these two cases. 3. Results and Discussion In this section, we present the results of the application of the Smoluchowski equation approach described in Section 2 to

a number of model systems. As described above, these model systems are based, in part, on Monte Carlo and molecular dynamics simulations of a model dye molecule dissolved in solvents (CH3I, CH3CN, and CH3OH) confined in spherical hydrophobic cavities of radius 10 and 15 Å. The simulations provide a nearly complete set of parameters for the Smoluchowski equation, albeit with some uncertainties, allowing a validation of the Smoluchowski equation approach. However, it is also interesting to examine the effect of the geometry of the confining framework (e.g., spherical vs cylindrical), the position dependence of the solvent reorganization force constant, k(r), and the position-dependent solute diffusion constant, D(r). We present the results of such an investigation here and discuss the implications. 3.1. Spherical Cavities. We first discuss results for timedependent fluorescence in hydrophobic spherical cavities, for which previous Monte Carlo and molecular dynamics simulations have been carried out to obtain free-energy surfaces and time-dependent Stokes shifts.17,35 The solution of the Smoluchowski equation gives the full two-degree-of-freedom timedependent distribution, f(r, ∆E, t). However, it is convenient to consider the two projections of the distribution, the distribution of solute radial positions, fr(r, t) (defined by eq 14) and the solvent coordinate distribution

f∆E(∆E, t) )

∫0R f(r, ∆E, t)dr

(16)

These two distributions are plotted for several times after excitation of a solute in CH3I confined in a 10 Å-radius cavity in Figure 5. It is also useful to examine the trajectory of the distribution mean, defined by

〈r(t)〉 )

∫0R dr ∫-∞∞ d∆E r f(r, ∆E, t)

(17)

∫0R dr ∫-∞∞ d∆E ∆E f(r, ∆E, t)

(18)

and

〈∆E(t)〉 )

This trajectory is plotted for the same system in Figure 6 in terms of 〈d(t)〉 ) R - 〈r(t)〉 and 〈∆E(t)〉. Examination of Figures 5 and 6 reveals a number of interesting features. First, we can note that the distribution of solute positions shown in Figure 5a has three peaks that correspond to, in decreasing distance from the cavity center, the solute lying next to and parallel with the cavity wall, the solute abutting but perpendicular to the cavity wall, and the solute roughly parallel to but away from the wall.35 In the ground state, the solute has a small dipole moment and sits primarily next to the wall in a parallel orientation, whereas in the excited state it is more evenly distributed with the free-energy minimum favoring positions near the cavity center but with more space available (∝4πr 2) near the spherical cavity wall. Second, it is clear that, as noted above, the solvent reorganization is rapid relative to the solute diffusion for these systems. Note that a shorter time interval is used for plotting f∆E(∆E, t) than fr(r, t) in Figure 5 to capture this faster motion. The effect is also evident from the trajectory shown in Figure 6 where the initial motion after excitation is primarily in the solvent coordinate; a significant fraction of the solvation response comes from this relaxation occurring in the first few picoseconds. Third, this solvent reorganization is followed by a slow solute diffusion away from the cavity wall that occurs on the time scale of tens to hundreds of picoseconds. From Figures 5a and 6, it can be seen that the solute does diffuse in the first 10-20 ps, but it is


Feng and Thompson

Figure 5. Time evolution of (a) radial distribution fr(r, t) and (b) the fluorescence spectrum, f∆E(∆E, t) of a solute in CH3I in a spherical cavity with R ) 10 Å, Case I.

Figure 6. Trajectory of the distribution mean of f(r, ∆E, t) is shown for a 10 Å-radius cavity containing CH3I; the average distance from the cavity wall, 〈d(t)〉 ) R - 〈r(t)〉, is plotted against the average solvent coordinate, 〈∆E(t)〉. Results for Cases I (black), II (red), III (blue), and IV (green) are shown. The time interval between points is 1 ps for the first 20 and 10 ps thereafter; the initial point is in the lower-right quadrant.

primarily diffusion from the parallel position next to the wall to the perpendicular orientation relative to the wall; the diffusion into the cavity interior takes place on a longer time scale. This diffusion is naturally accompanied by a simultaneous solvent reorganization, but with a magnitude smaller than that for the first several picoseconds. Fourth, a nonmonotonic change in the solvent coordinate can be seen for Case IV in Figure 6. This is because of the choice for k(r), which decreases with the distance from the cavity center (see Figure 1); this choice is likely unphysical for most real systems. Finally, it should be noted that the ultrafast inertial solvation, which lasts for a few hundred femtoseconds or less, predicted by molecular dynamics simulations and observed in ultrafast measurements of TDF, is not present in the results; the Smoluchowski equation is not applicable on that time scale.21 In the following, we discuss the characteristics of the solute diffusion after excitation and the time-dependent fluorescence. The effect of changing the cavity size and the solvent are specifically examined. 3.1.1. Solute Diffusion Dynamics. The diffusion dynamics of the solute can be characterized by calculating the average radial

Figure 7. Radial displacement of the solute, 〈∆d(t)〉, is plotted for CH3I solvent in a spherical R ) 10 Å cavity; the short-time behavior is plotted in the inset. Results are shown for Cases III (solid black line) and VI (dashed black line); 〈∆d(t)〉 for the other cases lies between these two curves. Molecular dynamics simulation results17 (solid red line) are plotted for comparison.

displacement of the solute away from the cavity wall,17 〈∆d(t)〉 ) 〈d(t)〉 - 〈d(0)〉. Results for 〈∆d(t)〉 obtained from solving the Smoluchowski equation for CH3I in a spherical cavity with a radius of 10 Å are shown in Figure 7 for Cases III and VI (the other cases lie in between the results for these two). The time scales of this radial diffusion can be extracted by fitting 〈∆d(t)〉 to a multiexponential function. Usually, 〈∆d(t)〉 can be fit reasonably with a biexponential function; however, a better fit is generally obtained with a three-exponential function

〈∆d(t)〉 ) Ad1[1 - e-t/τd1] + Ad2[1 - e-t/τd2] + Ad3[1 - e-t/τd3] (19) with three well-separated time scales. Here, the weights, Adj, denote the contributions to the radial displacement occurring with each time scale, τdj; the parameters are ordered in increasing magnitude of the time-scale, that is, τd3 > τd2 > τd1. The longest time scale, τd3, which is primarily associated with the radial solute diffusion from near the wall to the cavity interior, is not sensitive to the form of the fitting function; similar values are obtained from a biexponential fitting. It is interesting that for



TABLE 2: Values from Fitting the Radial Displacement 〈∆d(t)〉 to the Three-Exponential Fitting Function Defined by Equation 19 (Amplitudes in Å and Time Constants in ps) for the Solvent CH3I R (Å)

Case

Ad1

τd1

Ad2

τd2

Ad3

τd3

10 10 10 10 10 10 15 15 15 15 15 15

I II III IV V VI I II III IV V VI

0.011 0.013 0.015 0.020 0.012 0.015 0.039 0.032 0.054 0.012 0.039 0.033

3.8 3.5 2.6 4.3 4.1 4.8 4.4 4.9 3.9 6.6 4.1 3.8

0.097 0.099 0.10 0.091 0.10 0.10 0.15 0.16 0.14 0.13 0.26 0.35

14.6 14.5 14.5 14.5 13.7 13.7 46.8 45.0 43.9 41.4 39.4 39.1

1.39 1.40 1.39 1.40 1.39 1.38 2.14 2.11 2.18 2.09 2.02 1.94

156 156 157 155 153 150 70.4 70.0 70.4 71.9 69.1 69.8

TABLE 3: Values from Fitting the Radial Displacement 〈∆d(t)〉 to the Three-Exponential Fitting Function Defined by Equation 19 (Amplitudes in Å and Time Constants in ps) for the Solvent CH3CN R (Å)

Case

Ad1

τd1

Ad2

τd2

Ad3

τd3

10 10 10 10 10 10 15 15 15 15 15 15


0.021 0.022 0.023 0.029 0.020 0.020 0.038 0.032 0.052 0.014 0.039 0.030

3.7 3.6 3.0 4.7 3.6 3.7 4.6 5.2 3.8 7.8 4.3 3.7

0.097 0.096 0.097 0.090 0.10 0.11 0.87 0.88 0.89 0.80 1.02 1.10

22.5 22.5 22.5 22.8 20.6 20.4 54.1 53.8 53.9 52.5 50.7 51.0

1.13 1.13 1.13 1.13 1.12 1.12 1.87 1.86 1.88 1.90 1.73 1.65

228 228 229 228 222 215 136 136 137 147 138 142

TABLE 4: Values from Fitting the Radial Displacement 〈∆d(t)〉 to the Three-Exponential Fitting Function Defined by Equation 19 (Amplitudes in Å and Time Constants in ps) for the Solvent CH3OH R (Å)

Case

Ad1

τd1

Ad2

τd2

Ad3

τd3

10 10 10 10 10 10 15 15 15 15 15 15


-0.019 -0.021 -0.018 -0.027 -0.021 -0.021 0.13 0.12 0.15 0.084 0.13 0.11

3.4 3.4 4.5 3.7 2.9 2.9 6.4 6.9 6.0 8.6 5.9 5.6

0.61 0.61 0.64 0.56 0.63 0.64 0.48 0.49 0.48 0.47 0.65 0.78

19.9 19.8 20.0 19.4 17.7 17.5 44.7 44.9 43.9 43.1 42.2 42.5

1.75 1.75 1.76 1.72 1.72 1.71 2.06 2.03 2.08 2.01 1.88 1.77

139 138 140 137 136 132 89.8 89.7 89.8 95.0 89.1 90.6

the first few picoseconds, which involve primarily solvent reorganization, 〈∆d(t)〉 shows a slow nonexponential increase (see the inset in Figure 7). This is due to the decrease in the barrier to radial diffusion associated with the relaxation in ∆E. This nonexponential part contributes little to the radial displacement and is removed before fitting. The weights and time scales obtained from fitting 〈∆d(t)〉 are shown in Tables 2-4. It should be noted that the faster time scales, τd1 and τd2, have relatively small weights (Ad1 in particular) and thus these two time scales are not determined as accurately as τd3. The molecular dynamics simulation results plotted for comparison in Figure 7 are in good agreement with the Smoluchowski equation results; this is naturally expected because the diffusion coefficient was obtained by fitting to that data. However, the time constants reported in Table 2 of ref 17 by fitting these MD results are not as comparable: three time-scales are not resolvable and, while qualitatively accurate, the time scales and amplitudes are not

Figure 8. Plot of inverse of the longest time scale, 1/τd3, is plotted vs the diffusion constant, D. Results are shown for spherical cavities containing CH3I with R ) 10 Å (black line with circles) and R ) 15 Å (black line with squares), CH3CN with R ) 10 Å (red line with circles) and R ) 15 Å (red line with squares), and CH3OH with R ) 10 Å (blue line with circles) and R ) 15 Å (blue line with squares). The data reported in Tables 2-4 are indicated by an × .

as reliable when obtained from 100 ps dynamics rather than the 1 ns dynamics obtained with the Smoluchowski equation approach. The three time constants obtained from fitting 〈∆d(t)〉 are consistent with the discussion given above of the time evolution of fr(r, t), for example, as shown in Figure 5a. In particular, the three time scales can be understood by inspection of fr(r, t). The shortest time constant, τd1, is on the order of ∼2-8 ps and is associated with a small amplitude, Ad1 < 0.04 Å for CH3I and CH3CN and Ad1 e 0.15 Å for CH3OH (Ad1 is negative for CH3OH and R ) 10 Å, indicating the net effect is a rearrangement to positions closer to the wall). It is on this relatively fast time scale that the equilibrium ground-state distribution rearranges within each excited-state local minimum after excitation (i.e., it involves minimal barrier crossing), leading to only small changes in the average distance from the cavity wall. The intermediate time constant, τd2, is ∼13-20 ps for 10-Å-radius cavities and ∼40-55 ps for 15-Å cavities. This time scale corresponds primarily to diffusion over the smaller barrier separating parallel and perpendicular configurations of the solute adjacent to the wall. This can also be seen in Figure 5a and naturally corresponds to a larger amplitude with Ad2 typically at least three times Ad1. The longest time scale, τd3, is associated with the solute migration from positions near the cavity wall to the cavity interior and thus has the largest amplitude (Ad3 is between 1.1 and 1.8 Å for R ) 10 Å and 1.6-2.2 Å for R ) 15 Å). The time scale is due to both the longer migration distance and the barriers to solute diffusion. It is interesting that this three-stage picture of the diffusive process is qualitatively the same for the three solvents and two cavity sizes examined in this work. From Tables 2-4, it is clear that the longest time scale, τd3, decreases when the cavity radius is increased from 10 to 15 Å. This is in general agreement with experimental observations and previous nonequilibrium MD simulations.17 In the present Smoluchowski equation approach, we can identify this change with the increase in the solute diffusion coefficient. This is demonstrated in Figure 8 where 1/τd3 is plotted versus D for CH3I, CH3CN, and CH3OH in 10- and 15-Å cavities. Clearly, the two quantities are linearly related. For CH3I, when the value of D for the R ) 15 Å cavity (9.0 × 10-6 cm2/s) is used in a simulation for the 10 Å cavity, the longest time scale decreases

18068 J. Phys. Chem. C, Vol. 111, No. 49, 2007 from 156 to 79 ps, which is close to the value of 70 ps found for the larger cavity. This indicates that for CH3I the decrease in the longest time scale with R comes mainly from the increase in the diffusion constant. The changes in the free-energy surface have a lesser effect on the longest time scale. However, for the other two solvents, CH3CN and CH3OH, the changes in the freeenergy surface play a more important role than those in CH3I, and contrary to the CH3I case, the changes in the free-energy surface increase the longest time scales. Note the qualitatively different results for CH3CN and CH3OH compared to CH3I. Although for the latter τd3 decreases upon going from R ) 10 to 15 Å for a fixed value of D, the former show the opposite trend. The relative independence of τd3 on the free-energy surface in r occurs despite the fact that the free-energy curves change significantly with cavity size. The number of barriers and their heights are dependent on R with the free-energy curves becoming flatter as R increases.35 This does have a large effect on the second-longest time scale, τd2. As shown in Tables 2-4, τd2 is more than doubled when the radius is increased from 10 to 15 Å. This is due to the changes in the free-energy curves. Specifically, because the free-energy curve becomes flatter, the diffusion occurring on the intermediate time scale involves a greater migration in a 15-Å cavity than in a 10-Å cavity. This leads to a longer time for completion of this diffusive rearrangement. Because solvent reorganization happens on the order of a few picoseconds (see below), it does not significantly affect the longtime decay components of 〈∆d(t)〉. Rather, as noted above, the differences between the time scales for different solvents come mainly from different diffusion coefficients and free-energy curves for solute diffusion. For CH3CN, the longest time scales are 228 and 116 ps in 10 Å and 15 Å cavities, respectively. If the diffusion coefficients for CH3CN are used for CH3I, then the longest time scales for CH3I (156 and 70 ps for 10 Å and 15 Å cavities, respectively) become 266 and 109 ps. This indicates that for the 15 Å cavity the difference between the longest time scales for CH3I and CH3CN comes mainly from the difference in diffusion coefficients. In the case of the 10 Å cavity, both the changes in diffusion coefficients and in free energies contribute significantly. Note that there is a competition between the effects because of D and the free energy. That is, as the cavity is made larger, the distance the solute must diffuse is greater leading to slower time scales while the diffusion coefficient increases leading to shorter time scales. Because we do not have data on diffusion coefficients for CH3OH in these nanocavities, we have used the same values as those for CH3CN. Therefore, the differences in diffusion between the two solvents come completely from the differences in free-energy curves. The long (τd2 and τd3) time scales drop by about 60% from CH3CN to CH3OH. For 10-Å spherical cavities, the weight associated with the longest time scale for CH3OH is much larger than that of CH3CN because the difference between the initial and equilibrium distributions is much larger than that of CH3CN. The excited-state free-energy barrier around the highest peak of initial distribution for CH3OH is smaller than that for CH3CN; this would make the intermediate time scale faster. However, this diffusion involves a longer migration in CH3OH, so τd2 is actually similar for the two solvents. 3.1.2. Time-Dependent Fluorescence. The time scales for solvation dynamics are frequently extracted from measurements of the time-dependent Stokes shift. The result is generally plotted in the normalized form

Feng and Thompson

S(t) )

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

(20)

where 〈∆E(t)〉 is the fluorescence energy at time t after excitation. In the Smoluchowski equation approach, it can be obtained according to eq 18. The S(t) obtained from solving the Smoluchowski equation for CH3I in a spherical cavity with R ) 10 Å is shown in Figure 9a. From Figure 9a, one can see that the solvation comes mainly from the fast solvent reorganization, so the time scale related to this process can be obtained accurately. The slowest component can also be obtained reliably if the contribution to solvation from this process is not too small. The other time scales have small contributions to the normalized Stokes shift and thus may be subject to large errors. We have found that it is often sufficient to fit the normalized Stokes shift to a three-exponential function

S(t) ) Ae1 e-t/τe1 + Ae2 e-t/τe2 + (1 - Ae1 - Ae2) e-t/τe3

(21)

or a four-exponential function

S(t) ) Ae1 e-t/τe1 + Ae2 e-t/τe2 + Ae3 e-t/τe3 + (1 - Ae1 - Ae2 - Ae3)e-t/τe4 (22) The three-exponential fitting is used whenever two close time scales (differ by a factor less than two) are found in the fourexponential fitting. The time scales are arranged in such a way that τe1 < τe2 < τe3 < τe4 and each coefficient represents the contributions of the dynamics occurring on the corresponding time scale to the Stokes shift. The longest time scales cannot be extracted accurately by fitting S(t) to biexponential or threeexponential functions if its contribution to solvation is very small. Occasionally, a five-exponential function is needed for extracting the longest time scale. The fit parameters are given in Tables 5-7. From these multiexponential fit parameters, it can be seen that although some of the time scales are the same as those found for the diffusive motion represented by 〈∆d(t)〉, others are not. For example, the longest time scale of the Stokes shift decay is very close to the longest one for the diffusive motion. The second-longest time scale of decay of S(t) is on the same order of, but slightly smaller than, τd2, and both have relatively small amplitudes. Alternatively, the shortest time scale, τe1, is ∼1-4 ps, which is faster than the shortest time scale for diffusive motion (see Tables 2-4) and thus corresponds to solvent reorganization. This is also clear from Figure 6, which shows that the initial, rapid changes after excitation involve movement in the solvent coordinate, ∆E, with little change in the radial position. These results indicate that for the smooth hydrophobic nanocavities considered here the slow decay components of the Stokes shift are determined primarily by the diffusion of the solute toward the cavity center; the multiple barriers crossed in the diffusion lead to multiple-exponential decay in S(t). Because the time scale of solvent reorganization is fast, the solvent reorganization can follow the solute motion nearly adiabatically when the solute is diffusing toward the cavity center. Comparing Cases I-VI, the amplitudes of the long-time decay components can be understood based on the change in the position of the free-energy minimum in ∆E with the radial position r in the excited state as shown in Figure 2. Specifically, the degree to which the fluorescence energy is affected by the solute position can be seen from the angle the line connecting the minima in ∆E makes with the vertical axis in Figure 2. The larger this



Figure 9. Normalized Stokes shift, S(t), of a solute in CH3I in a spherical cavity with radius equal to 10 Å (left panel); the results, except Case IV, are plotted on a log scale in the inset. The Case I, V, and VI S(t) values are overlapping. The right panel shows the comparison of S(t) from a 200 ps MD simulation on CH3CN in a 15 Å cavity with the results from the Smoluchowski equation; the first 100 ps results are plotted on a log scale in the inset.

TABLE 5: Values from Fitting the Normalized Time Dependent Stokes Shift, S(t), Using Equations 21 or 22 for the CH3I Solvent R (Å) Case 10 10 10 10 10 10 15 15 15 15 15 15


Ae1

τe1

Ae2 (ps)

0.89 0.92 0.78 1.02 0.89 0.89 0.74 0.89 0.56 0.83 0.74 0.74

1.7 1.8 1.1 3.6 1.7 1.7 1.2 1.8 0.55 3.1 1.2 1.2

0.02 0.01 0.03 0.01 0.02 0.02 0.06 0.02 0.09 0.08 0.06 0.06

τe2

Ae3

τe3

12.5 0.09 155 12.9 0.07 154 10.7 0.19 156 11.0 -0.03 154 12.0 0.09 152 12.0 0.09 149 3.5 0.05 31.6 7.2 0.04 35.1 2.0 0.04 23.5 6.5 0.07 35.6 3.4 0.06 32.0 3.4 0.07 32.8

Ae4

τe4

0.15 0.05 0.31 0.02 0.14 0.13

69.0 69.6 68.7 71.6 67.8 68.4

TABLE 6: Same as Table 5 but for the CH3CN Solvent R (Å) Case 10 10 10 10 10 10 15 15 15 15 15 15

I II III IV V VI I II III IV1 V VI

Ae1

τe1

Ae2 (ps)

0.91 0.94 0.82 1.02 0.91 0.91 0.72 0.86 0.54 0.79 0.72 0.71

1.8 1.8 1.1 3.6 1.8 1.8 1.2 1.8 0.57 3.1 1.2 1.2

0.02 0.01 0.03 0.005 0.02 0.02 0.03 0.02 0.05 0.07 0.04 0.04

τe2

Ae3

τe3

14.3 0.07 214 14.4 0.05 214 11.4 0.15 215 30.7 -0.025 213 14.2 0.07 210 14.5 0.07 205 6.3 0.14 50.2 10.7 0.075 50.9 4.5 0.17 46.9 5.1 0.11 54.4 6.0 0.10 47.5 5.8 0.10 47.7

Ae4

τe4

0.11 0.04 0.24 0.03 0.14 0.15

132 133 132 145 134 137

angle, the larger the weights of the longest-time decay components because diffusion then has a greater effect on the fluorescence energy. If the alignment of the minima is nearly parallel with the vertical axis, then the amplitudes will be quite small (and the longest time scale will be difficult to extract from S(t)). It is important to note that the small amplitude components observed here may be missed in molecular dynamics simulations or experimental measurements because of statistical errors. For Case IV, the minima are angled away from the vertical axis, leading to the negative values of S(t) for some solvents at intermediate to long times. To our knowledge, such behavior has not been observed in experiments, indicating that the choice

Figure 10. Spectrum width of f∆E(∆E, t) obtained by fitting to a Gaussian distribution, 1/σx2π exp(-(((x - µ)2/2σ2))). The spectrum width is reported as 2.355σ. The results for Cases V and VI are nearly the same as those for Case I. The solvent is CH3I; the cavity is spherical with a radius equal to 10 Å.

TABLE 7: Same as Table 5 but for the CH3OH Solvent R (Å) Case 10 10 10 10 10 10 15 15 15 15 15 15


Ae1

τe1

Ae2 (ps)

τe2

Ae3

τe3

0.67 0.72 0.52 0.89 0.66 0.65 0.67 0.82 0.50 0.78 0.67 0.67

1.7 1.8 1.0 3.6 1.6 1.6 1.2 1.8 0.6 3.0 1.2 1.2

-0.03 -0.03 0.23 0.10 0.19 0.19 0.07 0.05 0.09 0.10 0.14 0.15

7.6 6.4 20.8 20.5 18.5 18.2 7.3 9.5 6.1 8.9 6.7 6.8

0.20 0.17 0.25 0.01 0.15 0.16 0.13 0.09 0.15 0.11 0.14 0.15

19.5 19.5 140 140 136 132 39.6 41.9 35.9 42.1 38.4 38.9

Ae4

τe4

0.16 138 0.14 138

0.13 0.04 0.26 0.01 0.11 0.10

87 88 87 88 87 88

of the force constant, k(r), for Case IV is not qualitatively realistic. In addition, previous Monte Carlo simulations for the systems considered here found that the width of fluorescence spectrum for the solute near the cavity center is wider than for


Feng and Thompson

Figure 11. Two slowest time scales extracted from the radial displacement, 〈∆d(t)〉 (left panel) and the normalized Stokes shift, S(t) (right panel). Results shown are taken from Case I. Each column represents the results for a particular solvent and geometry. The filled symbols denote results of 10 Å cavities and pores, and the open ones denote results of 15 Å cavities and pores. The circles represent results for the second-slowest time scales, the squares for the slowest time scales.

that near the cavity wall,18 meaning that k(r) decreases with decreasing r, consistent with Cases I and III and in contrast to Case IV. Comparing the results for Cases I, V, and VI provides information on the effect of the position dependence of the diffusion constant because these cases differ only in the parameter D(r). The results in Figure 9a and Tables 5-7 show that the position dependence of D(r) does not have a significant effect on the time-dependent Stokes shift. It is possible that such effects will appear for systems with stronger solvent-wall and/ or solute-wall interactions. To justify the Smoluchowski equation approach, we compared the normalized time-dependent Stokes shift from a 200-ps MD simulation on a solute in CH3CN in a 15 Å cavity with the results from the Smoluchowski equation in Figure 9b. As expected, the MD result shows an initial ultrafast drop from inertial motion of the solvent, which could not be observed in the Smoluchowski approach. After this initial drop, the MD result is in good agreement with the results from the Smoluchowski approach. Specifically, the longest time scale extracted from the MD result is 64 ps, and the one from Case II is 116 ps. The longest time scales usually have small amplitudes, which makes it difficult to extract long time scales from MD simulations. However, it is clear from the inset of Figure 9b that the long-time decay of S(t) predicted by the Smoluchowski equation is consistent with that obtained in the nonequilibrium MD simulations. The differences in the largest time constants presumably result from using a single time constant to fit the long-time component in the MD data but two for the Smoluchowski result as well as the shorter total time and statistical noise of the MD simulations. Interestingly, the amplitudes for the long-time component associated with the solute diffusion agree reasonably well between the MD and Smoluchowski results. This suggests that the large inertial component present in the MD simulations shows up as part of the ∼1-2 ps solvent reorganization component in the Smoluchowski equation approach, which cannot describe the inertial behavior. Finally, it is interesting to examine the time-dependent width of the fluorescence spectrum. From the evolution of the spectrum, f∆E(∆E, t), shown in Figure 5 it is apparent that the spectrum can be approximated well by a Gaussian distribution with a time-dependent center and width. An example of the evolution of the spectral width obtained from fitting the spectra

in this way is shown in Figure 10 for CH3I solvent with R ) 10 Å. From these results, we see that the spectral width increases monotonically with time, except for case IV, and the magnitude of the change is smallest for k(r) ) constant (Case II). The change in width can be attributed primarily to the difference between the ground-state and excited-state free-energy surfaces shown in Figure 2 because these free-energy surfaces determine the spectral widths of the initial and equilibrium distributions. For bulk solvents, the initial and equilibrium distributions have the same spectral width because the ground- and excited-state force constants are equal and not position-dependent. Carter and Hynes predicted an initial broadening of the spectrum followed by a narrowing in bulk solvents.32 This behavior was predicted to occur in less than 1 ps for their model system. Such nonmonotonic behavior at short time is not observed in the Smoluchowski equation approach; it means that this behavior is either system-specific or it is not captured by the Smoluchowski equation approach because of its inability to capture such short time dynamics. 3.2. Cylindrical Pores. One advantage of the Smoluchowski equation approach is that different confining framework geometries can be readily examined. We have used this to examine the differences between spherical and cylindrical confining geometries in the time-dependent fluorescence. It is almost certainly true that the free energies and diffusion coefficients are different for solvents in cylindrical pores than for solvents in spherical cavities. However, in the results presented here, the same free energies and diffusion coefficients are used for both frameworks to isolate the geometric effects. Thus, any differences between results for spherical cavities and cylindrical pores in this work come only from purely geometric factors. Our calculations indicate that the differences between spherical cavities and cylindrical pores due to only geometric factors are relatively small. An illustration is provided in Figure 11, which shows the time constant for the two slowest decay components for the radial displacement, 〈∆d(t)〉, and the normalized Stokes shift. Clearly, the geometry does not affect the time scales significantly. This is consistent with the results for spherical cavities for which the time scales for the longtime diffusive motion were found to be determined by the diffusion coefficient (see Figure 8) and the diffusive barriers. The largest difference is observed for the CH3CN solvent with


J. Phys. Chem. C, Vol. 111, No. 49, 2007 18071 pores are dominated by changes in the diffusion coefficients and free-energy surfaces (which include the effect of the solvent reorganization energy). However, we cannot exclude the possibility that the geometric effects might be somewhat larger for qualitatively different free-energy surfaces.

Figure 12. Total radial displacement of the solute in three solvents in cavities of different sizes and geometries. Each column represents results for a particular solvent and geometry. Circles and squares denote results for 10 and 15 Å cavities, respectively.

R ) 15 Å, where the cylindrical pore has time constants that are ∼30% longer than those for the spherical cavity. The changes in shorter time scales and amplitudes extracted from 〈∆d(t)〉 include (1) a decrease in Ad1 and Ad2 of up to 50% for most cases, although these amplitudes are small and may subject to large errors, (2) a significant increase in Ad3, with the largest increase being ∼0.6 Å for CH3CN in a 15-Å pore, and (3) in most cases the shortest time scales, τd1, change little, the largest change is a 2.1-ps increase from 5.2 ps for CH3CN in a 15 Å pore. Considering instead the normalized Stokes shift, the changes found for the amplitudes and shorter time scales are (1) small changes ( 5 Å, respectively; 3.54 ( 0.03 × 10-6 (cm2/s) and 4.54 ( 0.05 × 10-6 (cm2/s) are found for 0 < r < 3.6 Å and 3.6 < r < 6.4 Å, respectively. (41) Larger variations in the diffusion coefficient with dye molecule position may be expected in other systems with stronger interactions with the confining framework surface.

Smoluchowski Equation Description of Solute Diffusion Dynamics and

Recommend Documents