Solvation Dynamics in Acetonitrile: A Study Incorporating Solute

of Chemistry, Colorado State University, Fort Collins, Colorado 80523, U.S.A. ..... Wiley Interdisciplinary Reviews: Computational Molecular Scien...
0 downloads 0 Views 320KB Size
J. Phys. Chem. B 2005, 109, 3553-3564

3553

Solvation Dynamics in Acetonitrile: A Study Incorporating Solute Electronic Response and Nuclear Relaxation Francesca Ingrosso,† Branka M. Ladanyi,*,‡ Benedetta Mennucci,† M. Dolores Elola,‡ and Jacopo Tomasi† Dipartimento di Chimica e Chimica Industriale, UniVersita` di Pisa, Via Risorgimento 35, 56126 Pisa, Italy, and Department of Chemistry, Colorado State UniVersity, Fort Collins, Colorado 80523, U.S.A. ReceiVed: September 28, 2004; In Final Form: NoVember 26, 2004

The solvent reorganization process after electronic excitation of a polar solute in a polar solvent such as acetonitrile is related mainly to the time evolution of the solute-solvent electrostatic interaction. Modern laser-based techniques have sufficient time resolution to follow this decay in real time, providing information to be confirmed and interpreted by theories and models. We present here a study aimed at the investigation of the different steps involved in the process taking place after a vertical S0 f S1 excitation of a large size chromophore, coumarin 153 (C153), in acetonitrile, from both the solute and the solvent points of view. To do this, we use accurate quantum mechanical calculations for the solute properties within the polarizable continuum model (PCM) and classical molecular dynamics (MD) simulations, both equilibrium and nonequilibrium, for C153 in the presence of the solvent. The geometry of the solute is allowed to change in order to study the role of internal motions in the time-dependent solvation process. The solvent response function has been obtained from the simulation data and compared to experiment, while the comparison between equilibrium and nonequilibrium MD results for the solvation response confirms the validity of the linear response approximation in the C153-acetonitrile system. The MD trajectories have also been used to monitor the structure of the solvation shell and to determine its change in response to the change in the solute partial charges.

1. Introduction The mechanism of solvation dynamics as a response to a change in solute-solvent interactions has received wide attention in the past three decades, both from experimental and theoretical points of view.1-7 Understanding the dynamics of solvation in polar solvents has important applications to reaction dynamics, especially for reactions involving electron and ion transfer. Since the late 1980s, theoretical and simulation studies have been peformed by many groups, as described in the review articles1-4,6,7 and references therein. The initial attention on solvent response and on the mechanism of solvation, conducted by using simple solute models, expanded to incorporate more realistic solutes and focused on modeling solute-solvent interactions. These developments have been motivated in part by improvement in the time resolution of experimental techniques that has allowed measurement of the ultrafast nondiffusive solvent response,5,6,8-10 which is more sensitive than the longer-time diffusive response to solute-solvent interactions. This work represents a further contribution to the study and modeling of solvation dynamics at a molecular level through a realistic, molecular description of both solute (coumarin 153, C153) and solvent (acetonitrile) by using a combination of quantum calculations and molecular dynamics (MD) methods. Acetonitrile has often been chosen as a prototypical smallmolecule aprotic polar solvent. Solvation dynamics in this * To whom correspondence should be addressed. E-mail: lamar.colostate.edu. † Universita ` di Pisa. ‡ Colorado State University.

bl@

solvent and its effects on solute electronic spectra have been extensively studied using different experimental techniques, such as time-dependent fluorescence8,11-20 and photon echo21-23 experiments. The time resolution of these experiments is now under 100 fs, short enough to detect the nondiffusive ultrafast response, which accounts for most of the decay in the solutesolvent energy gap.8,9,13,14 Theoretical studies of solvatochromism in acetonitrile have been carried out by means of computer simulations.24-26 Theoretical models,27-33 MD simulations34-38 and combined approaches39-42 have been applied to the dynamics of the solvation process. Solvation shell structure, electrostatic solutesolvent interactions, and the effect of the solvent on the electronic spectra of naphthalene43 and C15344 in acetonitrile have been investigated by combining semiempirical and classical MD approaches. Because of a large change in its dipole moment and a small change in its intramolecular structure on electronic excitation, C153 has been the chromophore of choice in many experimental studies of solvation dynamics.13,14,45-47 It is therefore worthwhile to improve the modeling of this chromophore in polar solvents and to implement these improvements into a solvation dynamics simulation. In previous work, the C153 geometries and partial charges to be used as parameters for the simulation were obtained from semiempirical or ab initio calculations in a vacuum.35,44,48 In this work, we perform quantum mechanical (QM) calculations for C153 in the presence of the solvent, by means of integral equation formalism (IEF)49 polarizable continuum model (PCM).50 The geometries for both the ground and excited states of C153 are optimized in acetonitrile, and

10.1021/jp0456032 CCC: $30.25 © 2005 American Chemical Society Published on Web 02/05/2005

3554 J. Phys. Chem. B, Vol. 109, No. 8, 2005

Ingrosso et al.

Figure 1. (a) Coumarin C153 structure and numbers used in this work to label atoms, and (b) PCM cavity.

the electrostatic potential (ESP) fit of the solvated system is used to obtain the partial charges. This allows us to achieve a better representation of the coumarin properties in solution. In addition, to further improve the solvation model in our simulation, we have included solute geometry relaxation. Even though the structure of C153 is considered rigid enough to have little effect on solvent reorganization, intramolecular changes have been hypothesized18,51 on the basis of the evaluation of the excited-state geometry. We have therefore included some of the internal degrees of freedom of the solute molecule, which in addition is allowed to translate and rotate in the simulation box. To be able to isolate the effects of the C153 flexibility, we have also run simulations in which we have constrained the molecular geometry, and we anticipate here that we will find differences in the solvation response in the longer time scale. The outline of this paper is as follows. Section 2 is devoted to a theoretical introduction to our work: in subsection 2.1, we describe the quantum mechanical method used for calculation of the C153 electronic structure in acetonitrile; subsection 2.2 contains the statistical mechanics background for the treatment of the simulation data, while in subsection 2.3 we discuss briefly some useful properties of the solvation shell. In section 3, we present and discuss the results that we have obtained from the QM calculations (section 3.1) and from the simulation (sections 3.2 and 3.3). The main conclusions for this work are drawn in section 4. 2. Theoretical Background In this section, we outline the most important theoretical procedures that we have used in our work. It is beyond the scope of this paper to give a complete and detailed exposition of the theories and their implementation, for which we refer the readers to the original works. 2.1. Ground- and Excited-State IEFPCM Calculations. Coumarins, and in particular C153 (Figure 1a), have been shown to have almost planar geometry and highly emissive intramolecular charge transfer in polar solvents.52 In addition, the transition from the ground state to the first excited singlet state in this molecule does not interfere with nearby transitions, giving rise to a charge redistribution which can be considered dipolar.13 Many previous theoretical works have focused on coumarins: quantum and semiempirical calculations of geometries and electronic properties have been performed35.44,51,53-55 in order to achieve a better understanding of solute properties related to electrostatic interaction with the solvent. In particular, C153 has larger masses and moments of inertia with respect to those of standard solvents. This property makes it an ideal probe for solvation dynamics, because these characteristics result in

a small contribution of solute motion to the solvation response.35,42 In the present work, we have used the IEF version49 of the PCM to study ground and excited states of C153 in acetonitrile. In the IEFPCM approach, the solvent is represented as a dielectric medium, characterized by its dielectric constant, and the solute is described as a QM charge distribution embedded in a molecular-shaped cavity within the dielectric (see Figure 1b for the graphical representation of this cavity for C153). The presence of the solute charge distribution polarizes the medium and gives rise to a solvent reaction field, described in terms of apparent charges located on the cavity surface. This reaction field is introduced in the Schro¨dinger equation determining the solute wavefunction in the presence of a polar solvent through a proper electrostatic operator, still defined in terms of the apparent charges: the details about the implementation and the solution of the resulting set of equations can be found in ref 49. The IEFPCM has been extended to treat the formation and relaxation of electronic excited states,56 and this extension has been exploited here. The exact description of all possible aspects involved in the time-dependent evolution of a solvated system after an electronic excitation is too complicated to model, but by introducing proper approximations, we can still achieve a reliable picture of solvation and obtain good results for the calculation of the properties of interest. First of all, we consider only solutesolvent electrostatic interactions. This is not a strong approximation in this case, due to the fact that both the solute and the solvent molecules are highly polar, and electronic excitation leads to a large change in the dipole moment of the C153, which is the main source of solvatochromic shift in solution. The additional approximation is to partition the solvation response to the perturbation induced by the electronic excitation in the solute into two terms. This partition is made on the basis of the consideration that, to a good approximation, the solvation response occurs on two distinct time scales. There is a fast portion of the response which is determined by the electrons of the solvent molecules: this component will instantaneously equilibrate with the electronic density of the solute in the excited state. By contrast, the part of the response related to the solvent nuclear and molecular motions, which is much slower, will initially remain frozen in its original equilibrium with respect to the ground state. This description corresponds to what is usually referred to as a “vertical nonequilibrium excitation” and to what we have called step 1 of the time evolution of the excited state. This initial situation of a partial solvent response (usually indicated

Solvation Dynamics in Acetonitrile

J. Phys. Chem. B, Vol. 109, No. 8, 2005 3555

as nonequilibrium response) will then evolve toward a complete equilibrium. This evolution can be modeled in different ways and is represented here as occurring in three steps. The subsequent two steps (step 2 and 3) are model situations, without a direct physical meaning, in which the solvent first fully equilibrates to the excited state electron density of the solute (step 2), and then, the solute relaxes to the geometry corresponding to the minimum of the excited state potential (step 3). In step 3, as in step 1, only the electronic part of the solvent response will change. Finally, in step 4, the final solvent equilibration to the solute excited state will be reached. In the IEFPCM framework, the different steps will correspond to calculations in which the solvent apparent charges are computed in terms of the optical (step 1 and step 3) or the static (step 2 and step 4) dielectric constant. By means of the IEFPCM calculations, we thus have the opportunity to describe some steps of the formation and relaxation of the excited state, but we are not able to continuosly follow the time evolution. We therefore use the QM calculation as a starting point for the MD simulation and a reference for a posteriori comparison with the simulation results. 2.2. Solvation Dynamics in Polar Solvents. An interesting aspect of polar solvation dynamics is the ultrafast regime of the solvation response, when the solvent molecules are small (water, acetonitrile). It has been shown through experiments10,12,13 and computer simulations10,34,35,42,57 that the most important portion of the decay of the solvation response in smallmolecule solvents such as water and acetonitrile occurs in a sub-picosecond time scale and that it is not generally possible to describe the functional shape of the decay in terms of a combination of exponential functions. In greater detail, the reorganization of the solvent molecules around the solute occurs through an initial inertial behavior, responsible for the initial Gaussian shape of the decay, which then propagates through the solvent structure by librational motion, eventually becoming diffusive. The diffusive tail can be single-exponential or a combination of exponential functions, depending on how many different diffusive relaxation modes are active during the longtime-scale evolution of the polarization. Polar solvation dynamics is most often measured by timeresolved electronic spectroscopy through the detection of the time-dependent Stokes shift of the fluorescence spectrum of a dye molecule.1,2,4,6 Fluorecence spectra are recorded at various delay times t after electronic excitation, which occurs at t ) 0. The solvation response is usually reported in terms of the time evolution of ν(t), the maximum in the fluorescence signal:

S(t) )

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

(1)

with ν(∞) corresponding to the maximum of the steady-state emission spectrum. Typically, the dye molecules are much larger than the solvent molecules and are structurally fairly rigid, so that S(t) evolves in time primarily as a result of solvent reorganization in response to the change in the solute electronic structure. The fluorescence spectrum shifts further to the red as time progresses. The total shift, ν(0) - ν(∞), measures the extent change in the solvation free energy due to the solute electronic transition. Large shifts are therefore observed for solutes such as C153 that undergo a large change in dipole moment in highly polar solvents.13 The use of MD simulations is aimed at reproducing the solvent response to the perturbation induced by the solute. The straightforward way to model this situation is to perform nonequilibrium simulations. Starting from independent equi-

librium configurations for the solvent in the presence of the ground-state solute, the solute charge distribution is changed to that corresponding to the excited state. Simulations are then run for a set of nonequilibrium trajectories, and the results are averaged over these trajectories. Under the assumption that the variation of the solute geometry in the excitation process is negligible (vertical transition), the fluorescence Stokes shift is due to the change, ∆E(t), in the solute-solvent potential34 that occurs as a result of the solute electronic excitation. The solvation response can be written as57

SRF(t) )

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

(2)

where the overbar denotes a nonequilibrium average in the perturbed system. 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. In this case, ∆E is electrostatic and can be expressed as57

∆E(t) )



∆qRVR

(3)

R∈solute

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

VR )

qjβ

∑ ∑ j)1 β∈solvent 4π r

(4)

0 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. If a perturbation of a system at 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.58 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 fluctuations, δ∆E ) ∆E - 〈∆E〉, in the equilibrated unperturbed system (i.e., the solvent in the presence of the ground-state solute)57

SRF(t) = C0(t) )

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

(5)

where 〈‚‚‚〉0 indicates that the equilibrium ensemble average is performed in the unperturbed system. The range of validity of the linear response for solvation dynamics can further be tested by comparing the nonequilibrium response function to the time correlation

C1(t) )

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

(6)

of δ∆E in the equilibrated system containing the excited state (S1) solute in the presence of the solvent, with the equilibrium average in this system denoted as 〈‚‚‚〉1. 2.3. Properties of the Solvation Shell. In an isotropic fluid, the molecular pair distribution function can be decomposed into a rotational invariant expansion in which the coefficients are related to orientational correlations in the fluid. In the present case, we focus on pair correlations that describe the relative orientations of vectors embedded in the solute and the solvent

3556 J. Phys. Chem. B, Vol. 109, No. 8, 2005

Ingrosso et al.

dipole. For orientational correlations of this type, the coefficients hmnl(r), where the three indices (m, n, l) are related to the dependence of the pair correlation on three unit vectors, two embedded in the molecules and one along the intermolecular distance,59 are given by

mnl

h

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

(7)

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. We focus our attention on rotational invariants Φ110(12) and Φ112(12), which are most closely related to intermolecular dipolar correlations:

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

(8)

Φ112(12) ) 3(uˆ 1·rˆ12)(uˆ 2·rˆ12) - uˆ 1·uˆ 2

(9)

Φ110(12) is related to the mutual orientation of the molecules, while Φ112(12) has the angular dependence of the dipole-dipole interaction between them. Both depend on unit vectors uˆ 1 and uˆ 2 lying along particular axes embedded in molecules 1 and 2: In our case, one is along the acetonitrile symmetry axis, and the other one along either the C3-C8 or N13-C2 direction within the solute (see Figure 1a and section 3.1). Φ112(12) also depends on the unit vector rˆ12 along the intermolecular distance between the solvent center-of-mass and the midpoint of a solute intramolecular vector. The reason for the choice of the solute C3-C8 and N13-C2 intramolecular directions is to test the hypothesis about charge transfer occurring along these two directions, on the basis of the QM calculations. We also calculate the normalized function

G1(r) )

h110(r) 3g000(r)

) 〈cos θ12〉

(10)

where the coefficient g000(r) ) h000(r) + 1 is the radial distribution connecting the solvent molecule center-of-mass to the midpoint of a solute intramolecular vector. G1(r) represents the average of Φ110(12) ) uˆ 1‚uˆ 2 ) cos θ12 at intermolecular separation r.60 The results obtained from the calculation of properties defined in eqs 8-10 are discussed in section 3.3. 3. Results 3.1. Quantum Mechanical Calculations. The C153 solute molecule is a relatively large system to be treated at a QM level. The large number of electrons to be considered for the evaluation of the electronic states involved in the transition imposes a limit on the level of calculation that can be performed, as well as on the quality of atomic basis set that can be used.35,51,53,54 As a start, we have performed a geometry optimization in the ground state, both in a vacuum and in acetonitrile. In both cases, we have employed a 6-311G basis set including the d functions at the density functional theory (DFT) level, using the nonlocal exchange correlation functional by Becke, Lee, Yang, and Parr61 (B3LYP). The geometry in a vacuum and in solution for the first electronic excited state has been optimized at the CIS62 level using the same basis set, which has also been used for all the calculations presented in the following.

TABLE 1: Radii for the Spheres Used in Construction of the Cavity for C153 atom/group

R′ (Å)

C (1, 3, 4, 5, 6, 8, 10, 33) C-H (2, 11) C-H2 (14, 15, 18, 23, 24, 27) N (13) O (9, 12) F (34, 35, 36)

1.7 1.9 2.0 1.55 1.52 1.47

TABLE 2: Excitation Energy and Ground-State Dipole Moment for C153 in Vacuum ∆ETD-DFT ) 3.24 eV quantitya

DFT

HF

µx (D) µy (D) µz (D) |µ| (D)

-3.41 -1.37 -6.75 7.68

-3.45 -1.21 -5.88 6.92

a Dipole moment components in this table and in Tables 3-5 are in debye units (1 D ) 3.3356 × 10-30 C m).

TABLE 3: Excitation Energy and Ground-State Dipole Moment for C153 in Acetonitrile ∆ETD-DFT ) 3.05 eV quantity

DFT

HF

µx (D) µy (D) µz (D) |µ| (D)

-4.81 -1.81 -9.08 10.44

-5.05 -1.67 -8.09 9.68

The calculations in acetonitrile have been performed by using the IEFPCM described in section 2.1. In particular, for both ground and excited states, we have used the cavity resulting from the overlap of spheres centered on atoms or groups of atoms. The radius of the spheres is equal to 1.2× the van der Waals atomic radius.63,64 In Table 1, we present the cavity details, while the cavity is shown in Figure 1b. All of the QM calculations were carried out using the Gaussian 03 software.65 From the geometry optimizations in a vacuum and in solution, it has been found that C153 is always in a syn conformation of the julolidyl rings, both in the ground and excited states. We note that C153 has two conformational isomers, characterized by a syn or anti arrangement with respect to the nitrogen atom.53,54 To determine the vertical transition energy, ∆ETD-DFT, a timedependent-DFT (TD-DFT) calculation has been run both in a vacuum and in solution. For details about the TD-DFT methods and their implementation in the IEFPCM scheme, we refer the interested reader to ref 66 and references therein. We have also performed calculations of the dipole moment µ of the ground state (Hartee-Fock (HF), DFT) and of the partial ESP charges on the atomic centers of C153. The results obtained in a vacuum are collected in Table 2, while the calculations in the presence of the solvent are collected in Table 3. As for the excited-state calculations, we have studied how the dipole moment of the molecule changes with the evolution of the excited state. In a vacuum, we have the vertical transition (step A) and the final equilibrium excited state (step B). The results of the CIS calculations are collected in Table 4 (vacuum). As said in section 2.1, four steps have been considered in acetonitrile starting from the vertical transition (step 1) to the final equilibration (step 4). The results are collected in Table 5, and they also include the fictitious steps 2 and 3. Because we have performed CIS calculations in the excited state, the ground-state charges have been obtained from HF

Solvation Dynamics in Acetonitrile

J. Phys. Chem. B, Vol. 109, No. 8, 2005 3557

TABLE 4: Excited-State Dipole Moment of C153 in Vacuum

TABLE 6: Partial Charges on the C153 Sites Significantly Involved in Charge Transfer

quantity

step A

step B

site

qvac,gr (a.u.)

qvac,exc (a.u.)a

qac,gr (a.u.)

qac,exc (a.u.)b

µx (D) µy (D) µz (D) |µ| (D) ∆|µ| (D)a ∆θ (°)b

-4.38 -2.00 -9.71 10.84 3.92 6.0

-4.35 -2.00 -9.59 10.72 3.80 5.9

C1 C2 C3 C8 O9 C10 O12 N13 C23 C27 C33 F34 F35 F36

-0.16 -0.26 -0.13 -0.08 -0.44 0.96 -0.64 -0.54 0.69 0.64 1.17 -0.37 -0.36 -0.33

-0.03 -0.39 0.11 -0.34 -0.44 0.86 -0.62 -0.46 0.61 0.75 1.36 -0.42 -0.42 -0.35

-0.21 -0.23 -0.18 -0.04 -0.48 1.04 -0.73 -0.49 0.60 0.67 1.21 -0.38 -0.37 -0.34

-0.16 -0.33 0.04 -0.28 -0.42 0.84 -0.73 -0.36 0.53 0.73 1.33 -0.43 -0.42 -0.36

a Change in the dipole moment magnitude between excited and ground states. b Change in the orientation of the ground-state dipole moment relative to the excited state (represented as the angle between the two different directions).

TABLE 5: Excited-State Dipole Moment of C153 in Acetonitrile quantity

step 1

step 2

step 3

step 4

µx (D) µy (D) µz (D) |µ| (D) ∆|µ| (D)a ∆θ (°)b

-6.28 -2.63 -13.00 14.67 4.99 6.1

-6.67 -2.90 -14.40 16.17 6.49 7.0

-6.65 -2.83 -13.66 15.45 5.77 5.9

-6.93 -3.04 -14.86 16.68 7.00 6.9

a Change in the dipole moment magnitude between excited and ground states. b Change in the orientation of the ground-state dipole moment relative to the excited state (represented as the angle between the two different directions).

Figure 2. Charge difference (a.u.) between the excited and ground states of C153 in acetonitrile in the four steps considered.

calculations in order to have a consistent comparison. In Figure 2, we report the difference between the excited- and groundstate charges for the most significant atomic centers of C153 in acetonitrile in the four steps described already, while in Table 6, we show the result of the calculation of the partial charges on these atoms in the ground state and in the equilibrium excited state, both in a vacuum and in acetonitrile. The dipole moment of the molecule µ in the ground and equilibrated excited states is shown in Figure 3a, where a visual representation based on different colors for different values of the partial charges in both states is provided. The direction of the µ vector follows the dipolar distribution resulting from the alternation of positive and negative charges going from C14 to O12. In Figure 3b, the molecular frame axes are shown (the y-axis is pointing inward) together with the dipole moment. After excitation, the orientation of dipole in the x-z plane changes slightly, as illustrated in Figure 3, as can be seen from Tables 4 and 5 in which we also report the angle between µ in

a

Step B. b Step 4.

the excited and ground states, and the magnitude of its out-ofplane component increases. A large increase of the dipole magnitude is found upon excitation, and this increase is much larger in acetonitrile than in a vacuum. This strong effect due to the presence of the solvent can also be seen by studying how individual Cartesian components of µ change. In going from the ground to the excited state in a vacuum, the increase in µx is 26%, in µy 65%, and in µz 63% of the ground-state value, while in acetonitrile we have, respectively, 37%, 82%, and 84%. The common feature between vacuum and solution is that the z- and y-directions are the most involved in the change. A possible explanation is that the charge transfer between the donor (amino group) and the acceptor (-CF3 group) is mainly a shift of electronic density along the z-axis mediated by the aromatic electrons, located along the y-direction. An analysis of the charge difference through the different steps considered allows us to make a hypothesis on the mechanism of the charge transfer in C153. Referring to Figure 2, we can recognize two main charge fluxes, one between C3 and C8 and the other one between C2 and N13. We have tested this preliminary consideration with the aid of MD simulations, and we anticipate that the results that we have obtained are in agreement with it (see section 3.3 for details and quantitative discussion). 3.2. MD Simulation Details. We have performed simulations of one C153 molecule in 255 acetonitrile molecules at T ) (298 ( 7) K in the microcanonical ensemble using the CCP5 software DL_POLY_2.67 The equations of motion were integrated using the Verlet leapfrog algorithm.68 The SHAKE algorithm69 was used for intramolecular structural constraints. The Ewald summation method with conducting boundary conditions68 was used for the long-range interactions. The simulations were run with cubic periodic boundary conditions and the Lennard-Jones (LJ) potential was cut off at one-half the box length. The simulation box length was chosen to reproduce the experimental density 0.777 g/cm3 at atmospheric pressure at the temperature considered.70 The solvent (CH3CN) has been treated with an all-atom representation as a rigid molecule using the parameters in ref 71, while we have included some internal degrees of freedom for the solute molecule, also represented at an atomistic level. In particular, we have considered the aromatic portion of the C153 molecule (i.e., the two central rings) and all the -CH2 groups as rigid units, while we have allowed stretching, bending, and torsions in the rest of the molecule.

3558 J. Phys. Chem. B, Vol. 109, No. 8, 2005

Ingrosso et al.

Figure 3. (a) Partial charges (a.u.) on atomic centers of the solute molecule in acetonitrile in the ground state and in the equilibrium excited state. The molecular dipole moment µ is represented as a red vector. The color scale used for the visualization is shown on the right-hand side. In (b), the reference system used in the calculations is displayed (the y-axis is pointing inward) together with the µ vector.

The intramolecular potentials used for these motions are in the form

U(rRβ) )

ks (r - r0)2 2 Rβ

(11)

for elongation of the bond between atoms R and β

U(θRβγ) )

kb (θ - θ0)2 2 Rβγ

(12)

for bending of the angle connecting atoms R, β, and γ, and a combination of three cosine functions for torsion, namely

U(φ) )

A2 A1 (1 + cos φ) + (1 + cos 2φ) + 2 2 A3 (1 + cos 3φ) (13) 2

The parameters in eqs 11-13 that we have used for the simulation are the OPLS values in the literature.72-74 The intermolecular potential between two atomic centers R and β on different molecules was represented as a sum of LJ and Coulombic potential terms:

uRβ ) 4 xRβ

[(

) (

σR + σ β 2r

12

-

)]

σR + σ β 2r

6

+

qRqβ 4π0r

(14)

where R and σR are the LJ well depth and diameter associated with the site R, qR is the partial charge on site R, and as eq 14 indicates, LJ interactions between unlike sites are constructed using Lorentz-Berthelot combining rules.68 The partial charges and the geometry of C153 in the ground state have been obtained from the IEFPCM calculation in solution, as described in section 3.1. For the equilibrium excited state, we have used the optimized excited-state geometry and the partial charges corresponding to step 4, while for the nonequilibrium simulation, we have used the ground-state geometry and the charges corresponding to step 1. We recall that the nonequilibrium simulation is the representation of what happens in the solute-solvent system immediately after the electronic excitation S0 f S1 (i.e., after a vertical transition as in step 1). For each equilibrium simulation (ground and excited states), an equilibration of 1 ns has preceded the data acquisition (1 ns). The time step used is 3 fs. The ground-state equilibrium trajectories were saved every 3 ps and used as initial configu-

ration for the nonequilibrium runs. The average in eq 2 has been calculated over 500 nonequilibrium trajectories. In addition to the simulations, described in the preceding paragraph, in which a partially flexible model solute was used, we ran MD simulations in which the C153 intramolecular structure was kept fixed. In this case, we used a FORTRAN code developed in our group. For the rigid-solute simulation, after equilibration runs of 500 ps, long (1.2 ns) equilibrated trajectories performed in the NVE ensemble at room temperature have been saved for C153 in the ground and excited states. Lennard-Jones forces were cut off at one-half the box length, and long-range Coulomb interactions were treated using the standard Ewald sum method with conducting boundary conditions. Cubic periodic boundary conditions were used in systems composed of one C153 and 255 solvent molecules, with a simulation box length that reproduces the experimental density at room temperature. The equations of motion were solved using the Verlet algorithm with a time step of 3 fs; all intramolecular distances were constrained with the SHAKE algorithm.69 The solvent has been modeled using the same all-atom representation mentioned already.71 The rigid C153 was modeled using the same partial charges, geometry, and LJ parameters as those employed in simulations in which stretching, bending, and torsion were allowed. 3.3. Solvation Shell Structure. The study of site-site solute-solvent radial distribution functions is aimed at comparing the arrangement of solvent molecules around the ground and excited states of the solutes. In particular, it is reasonable to expect that major changes would occur around the sites for which the partial change ∆qR is large. In Figure 4, we display some of the solute-solvent radial distributions gRβ(r)’s for which we have observed significant changes in going from the ground to the excited state of C153. Four of the displayed functions involve the acetonitrile methyl C atom, denoted as Cm, and two of the nitrogen site, N, which act, respectively, as positive and negative charge centers. The Cm site does not actually carry a positive charge but acts as a positive center, because it is surrounded by three positively charged H atoms. In the case of the solute, sites N13, O12, and F34 carry negative partial charges, while C10, C23, and C33 are positively charged (see Figures 1 and 3). We see that, in all cases, the excited state gRβ(r) corresponds to an enhancement in the solvation structure. We now consider the individual gRβ(r)’s, starting with the top panel, where gC10-Cm(r) and gO12-Cm(r) are displayed. The first peak of gC10-Cm(r) increases as one might expect, given that the positive charge on the C10 atom decreases in the excited

Solvation Dynamics in Acetonitrile

Figure 4. C153-acetonitrile site-site distribution functions: most significant changes between the ground and excited states.

state and the interaction with the Cm acetonitrile site becomes less repulsive. We see that the first peak of gO12-Cm(r) increases substantially in going from the ground to the excited state, even though the partial charge on the O12 atom does not change (see Figure 3 and Table 6). However, O12 experiences an induced effect due to its proximity to the C10 atom, for which the change is much larger, and to the fact that it is more accessible to the solvent than the C10 site. Our result for gO12-Cm(r) is in agreement with the behavior of this function calculated in ref 44, where the position of the first peak was found to be 3.18 Å (3.2 Å in our case) and where the increase in its height in going from the ground to excited state was found to be 0.63 (0.7 in our case). As for the second row in Figure 4, the behavior of gC33-Cm(r) seems at first glance to be in disagreement with expectations based on the change in the partial charge of the C33 site between ground and excited states. The Cm local density increases even though C33 becomes more positive. This can be interpreted as

J. Phys. Chem. B, Vol. 109, No. 8, 2005 3559 an effect of the three adjacent fluorine atoms, which become more negative. We also show in the figure gF34-Cm(r), which illustrates the increase in the Cm density around one of the F atoms. A similar effect is observed in the bottom panel of Figure 4, where gN13-N(r) and gC23-N(r) are displayed. An enhancement in the local structure is observed for both of them. In the excited state, the C23 positive site charge becomes slightly smaller, ∆qC23 ) -0.07 e, while the N13 negative charge decreases in magnitude, ∆qN13 ) 0.13 e, becoming less negative. The N local density augmentation on C23 is induced by its close proximity to the N13 atom. As already mentioned in section 2.3, it is possible to study the mutual orientation of solute and solvent dipoles using the projections of the solute-solvent pair correlation along rotational invariants Φ110 and Φ112 defined in eqs 8 and 9. On the basis of the values of partial charges in the ground and excited states of C153 (see Figure 2 and Table 6), the directions C3 f C8 and N13 f C2 are related to the flux of charge from the donor (amino) group to the acceptor (-CF3). Specifically, if we define the dipole moments of these vectors as pointing from the less positive to the more positive partial charge of each site pair, the ground-state dipoles for both C3 f C8 and N13 f C2 are in the direction of the arrows. In each case, the dipole moment change, ∆µ, resulting from electronic transition is in the direction opposite to the arrows. Given the locations of the sites within the C153 molecule, these directions of charge flow are in the prevailing direction of the dipole change for the whole molecule. We have therefore chosen these two directions to study the projections of solute-solvent pair correlations mutually along the orientation of dipoles (function h110(r)) and the dipolar interaction between them (function h112(r)). The results are shown in Figure 5. As can be seen from the figure, the changes in h110(r) and 112 h (r) in going from the ground to the excited state of C153 are similar for C3 f C8 and N13 f C2. In both cases, orientational correlations are stronger in the excited than in the ground state. In the ground state, h110(r) is negative near contact, corresponding to antiparallel dipole orientation, and then exhibits a positive peak at a somewhat larger r. Both features are more

Figure 5. Projection of rotational invariants: the orientational behavior of the solvent molecules with respect to the directions along which the charge transfer occurs in the solute are shown.

3560 J. Phys. Chem. B, Vol. 109, No. 8, 2005

Ingrosso et al.

Figure 6. Pair correlation functions which describe the dipolar orientational order of solvent molecules with respect to the directions along which the charge transfer occurs in the solute, as defined in eq 10.

pronounced for N13 f C2 for which the ground-state dipole is larger. In the excited state, the first peak of h110(r) is positive and much larger. Its sign is consistent with a predominantly parallel relative orientation of these intramolecular vectors and the solvent dipole. It can be explained by noting that ∆µ < 0 along both C3 f C8 and N13 f C2 vectors and that the overall solute dipole in the excited state is in a direction that is close to opposite of the directions of these vectors. The behavior of h112(r) is more complex, given that it depends on three relative orientations instead of just one. In the case of perfect antiparallel alignment of solute and solvent intramolecular vectors, rˆ12 near contact is perpendicular to uˆ 1 and uˆ 2, which would lead to h112(r) > 0 when h110(r) < 0. We see that, indeed, near contact, h112(r) and h110(r) have opposite signs. More generally, attractive dipole-dipole interaction leads to h112(r) > 0. In the case of both C3 f C8 and N13 f C2, h112(r) < 0 in the first solvation shell and in most of the distance range between contact and r = 8 Å in the excited state. This indicates again that the dipolar interaction of the solvents with these solute site pairs is favorable for the intramolecular vector pointing in the direction opposite to the arrows. In Figure 6, we display the function G1(r), defined in eq 10. These data do not contain as much information about the strength of intermolecular orientational correlations but can be straightforwardly related to the value of the angle θ12 in different intermolecular distance ranges. We see that, in going from the ground to the excited solute state, the changes in G1(r) for the two intramolecular vectors C3 f C8 and N13 f C2 are in the direction of decreasing angle θ12 (i.e., increasing 〈cos θ12(r)〉) relative to the solute dipole in the r range within the first solvation shell. For the C3 f C8 vector in the ground state, G1(r) at contact corresponds to around 101° and at the first peak of h110(r) to around 88°, while for the excited state, the θ12 value at the peak of h110(r) decreases to about 68°. For the N13 f C2 vector, the groundstate contact and peak values of θ12 correspond to around 128° and 86°, while the excited-state peak value decreases to around 77°. 3.4. Flexible Geometry Effects. Most of the low-to-moderate-frequency internal motions in the solute molecule are allowed in our simulation. Considering the internal degrees of freedom also provides a way to achieve a more accurate representation of time-dependent phenomena involved in solvation and to include some aspects of the polarizability of the molecule. Specifically, the aspects that are included involve fluctuations and static distortions in the molecular shape in response to the solvent environment.

Figure 7. Density distributions of the solute molecular dipole in three different stages of the excitation process reproduced in the simulation: ground state, vertical transition (nonequilibrium), and excited state.

The indirect inclusion of polarizability effects in the simulation can be shown with an example, by calculating the density distribution of the molecular dipole (Figure 7). For a rigid nonpolarizable solute molecule, the dipole in the molecular frame is constant throughout the simulation. In our case, these constant values, obtained from QM equilibrium geometries, correspond to the full vertical lines in the figure. It can be seen from the figure that the average dipoles for the flexible C153 are larger than for the rigid geometry, indicating that the molecule distorts to allow more favorable solutesolvent electrostatic interactions. The result that we obtain for flexible C153 is if not quantitatively (due to the different methods and model used) then qualitatively similar to the result obtained in ref 44, where the authors adopt a MD-semiempirical strategy to include the solute electronic polarization by the surrounding acetonitrile or methanol solvent. The C153 molecular dipole distribution that they obtain is peaked at a value larger than the gas-phase solute dipole. The dipole enhancement in that case is the result of the response of solute electron density to the fluctuations in local solvation structure. It would clearly be desirable to include the effects of both the solute shape and electronic polarizability fluctuations in a solvation model. The analysis that we have performed in order to determine if torsional motion contributes significantly to the dynamics is based on the investigation of the motion of the -CF3 group. One possible question to be asked is whether internal rotation occurs on the time scale of solvation dynamics. The dihedral C3-C8-C33-F angles have been calculated for each fluorine atom in the ground and excited (equilibrium,

Solvation Dynamics in Acetonitrile

J. Phys. Chem. B, Vol. 109, No. 8, 2005 3561

TABLE 7: Steady-State Stokes Shift: Comparison Between Calculations and Experiment source

Stokes shift (103 cm-1)

∆∆ν, eq 15 ∆∆νlra,0, eq 17 ∆∆νlra,1, eq 17 experiment, ref 13 Table 2 of ref 35 (∆∆νlra,0, eq 17)

2.0 ( 0.2 2.1 ( 0.2 2.2 ( 0.2 2.2 ( 0.2 3.2

nonequilibrium) states as functions of time. What we find is that there are only small oscillations around the equilibrium values of the torsional angles obtained from the QM geometry optimization, both in the ground and excited states. More information can be gained by studying the umbrella motion of the -CF3 group, which can be defined in terms of the improper dihedral angle between the plane defined by C33 and two of the fluorine atoms and the direction of the bond between C33 and the third F atom. Even in this case, the changes were only in terms of small displacements with respect to the equilibrium value of the angle. By combining the information derived from the two sets of angles, it is possible to conclude that the -CF3 group is not free to rotate on the time scale relevant to solvation dynamics. However, as will be shown in the next subsection, the effect of internal solute motion on the dynamics of solvation is nonnegligible and should be taken into account in order to achieve a more realistic description of the overall relaxation mechanism. 3.5. Steady-State Stokes Shift and Solvation Dynamics in Acetonitrile. The solvatochromic line shift of absorption and emission electronic spectra of a solute immersed in a polar solvent can be evaluated from MD simulations.35,38,75-77 Here, we used equilibrium simulations of both ground- and excitedstate C153 in acetonitrile. The steady-state Stokes shift ∆∆ν was calculated75,78 as the difference between the absorption and emission shifts:

∆∆ν ) ∆νabs - ∆νem

(15)

The shifts ∆νabs and ∆νem were estimated as averages of ∆E, obtained according to eq 3, for the ground- and excited-state equilibrium systems, respectively:

∆νabs(em) ) 〈∆E〉0(1)

(16)

A different way35 to calculate steady-state Stokes shifts from equilibrium simulation data is based on the assumption that the linear response approximation (lra) is accurate. Under this approximation, the line shift ∆∆νlra,0(1) is evaluated from the t ) 0 value of the ∆E fluctuation time-correlation function (see eqs 5-6) for the unperturbed system containing the solute in the ground (excited) state:

∆∆νlra,0(1) =

〈|δ∆E|2〉0(1) kBT

(17)

If the linear response approximation is valid, the two different calculations in the ground (∆∆νlra,0) and excited (∆∆νlra,1) states should give the same result. In Table 7, we compare the Stokes shift obtained in this work with the experimental result13 and with the result from another simulation of C153 in acetonitrile.35 Because the latter is the Stokes shift calculated according to eq 17 for the ground state, we also show the linear response results that we obtain for ∆∆νlra,0 and ∆∆νlra,1. These results are in agreement with the linear response prediction ∆∆νlra,0 ) ∆∆νlra,1 within the error bar associated with the calculation.

Figure 8. Solvation response function: the experimental function in ref 13 (dotted line) is compared to the result of the convolution (dashedand-dotted line) and to the fit as in eq 18 (dashed line).

TABLE 8: Calculated Parameters for Eq 18 Compared with the Experimental Ones13 calcd exptl (ref 13)

τ1 (ps)

% exp l

τ2 (ps)

% exp 2

0.09 0.09

82% 69%

0.86 0.63

18% 31%

The Stokes shift obtained in this work is smaller with respect to the one in ref 35 and in better agreement with the experimental result within the relative error bars. The dipole moment difference between excited and ground states that we have used is 0.9 D smaller than in ref 35, where the authors point out that scaling down the partial charges on the solute by 20% would lead to a good improvement in the Stokes shift result. It is worth mentioning that shifts of absorption and emission spectra for C153 in acetonitrile have also been evaluated44 using a combined semiempirical/MD methodology, which takes into account the solute polarizability. However, the procedure used to calculate the emission shift in ref 44 was different from the one used for the absorption shift. We have therefore not included the total ∆∆ν resulting from this calculation in Table 7. We finally present the results that we have obtained for the calculation of the solvation response. Its decay in acetonitrile is quite fast,13,35 because of the small dimensions of the molecule and its relatively small moments of inertia, as well as the highly collective character of the solvation response.42 To make a coherent comparison between the solvation response calculated from MD simulation and the experimental result, we have convoluted the former with a Gaussian function having the fullwidth at half-maximum equal to the experimental resolution13 (100 fs). Starting at t ) 100 fs, we have then performed a fit of the convoluted function to a sum of two exponentials

F(t) ) a1 e-t/τ1 + a2 e-t/τ2

(18)

the functional form used to fit the experimental data. In Figure 8, we show the result that we have obtained compared to the experimental function, on which the error bars of the experimental response function35 are displayed. Within the experimental error, the calculated and experimental functions are in extremely good agreement. In Table 8, we collect the numerical results of the fit procedure. It should be noted, however, that the short-time portion of the response obtained from MD is not well-represented by the functional form of eq 18, given that the inertial and librational portions of SRF(t) do not have the characteristics of exponential decay.

3562 J. Phys. Chem. B, Vol. 109, No. 8, 2005

Figure 9. Solvation response function: the calculated results (equilibrium ground state, dashed-and-dotted line; equilibrium excited state, dashed line; nonequilibrium, solid line) are shown.

The next step is to verify the validity of the linear response approximation. To do this, we have compared the solvation time correlations obtained from equilibrium simulations in the ground and excited states, given by eqs 5 and 6, to the SRF(t) obtained from nonequilibrium simulations. As shown in Figure 9, the three calculations of the solvation response give the same result, confirming that the lra accurately describes solvation dynamics in acetonitrile. This observation has already been made for simpler solutes34 as well as for a rigid model of C15335 in this solvent. Let us briefly discuss the validity of the lra in the system under study. Acetonitrile is a highly polar solvent with strong electrostatic interactions among solvent molecules. The perturbation represented by ∆E is not large compared to these interactions, even though it does lead to detectable restructuring of the solvation shell. This leads to good agreement with the linear response approximation. More generally, the lra has been found to work well for solvation dynamics in one-component polar liquid solvents7 as long as the solute-induced perturbation does not significantly change the pattern of solute-solvent hydrogen bonds35,79,80 or the solute size.81 The calculated functions feature a Gaussian shape at short times and librational oscillations that are not detected in the experimental S(t). It has already been pointed out35 that this is most likely due to the experimental resolution (0.1 ps) of the upconversion laser technique used to measure the time evolution of the fluorescence signal (see eq 1). To investigate the effect of considering a flexible coumarin in the MD simulation and to quantify the effect on its internal degrees of freedom on the solvation response, we have run simulations with C153 in a geometry that is fully constrained but free to translate and rotate in the simulation box. We compare in Figure 10 the results obtained for these two different C153 models. From the comparison, it is possible to notice that the inclusion of internal modes has an effect only on the longer time scales. The short-time dynamics is therefore dominated by acetonitrile motion, while the contribution due to the solute affects solvation dynamics on longer times. Solvation response in the rigid coumarin case is slower than in the flexible case. This suggests that solute flexibility provides additional relaxation channels for ∆E. We have applied the same procedure described here to compare the result for rigid C153 to the experimental result, and we show in Figure 11 the comparison, as in Figure 8.

Ingrosso et al.

Figure 10. Comparison between the SRF from two different simulations of C153 in acetonitrile: flexible coumarin (solid line) and rigid coumarin (dashed-and-dotted line).

Figure 11. As in Figure 8, the result of the convolution of the SRF obtained with rigid C153 in acetonitrile and the fit of the convoluted function (see eq 18).

As can be concluded by comparing the two results, the inclusion of internal motions in the solute appears to give an improved agreement with the experimental measurement. Summary and Conclusions In this paper, we have presented an investigation of timedependent solvation phenomena by means of a combination of different approaches. QM calculations of the solute electronic structure (C153) have been performed in the presence of the solvent (acetonitrile) and used both to study the S0 f S1 transition and to provide the parameters that we have employed for the MD simulation of solvation dynamics. We have run nonequilibrium and equilibrium (ground state, excited state) simulations of atomistic C153 in acetonitrile, including the potentials to describe solute intramolecular motions. The results of molecular dynamics have been compared with experimental data in the literature and discussed on the basis of our QM calculations. The results that we have obtained can be summarized in the following way: (i) The charge transfer between the donor and the acceptor group in the solute during the electronic excitation occurs along two main directions. With reference to Figure 1a, these two directions are identified by the vectors C3 f C8 and N13 f C2. As a consequence of the charge transfer, the properties of the solvation shell change. The changes concern the number of

Solvation Dynamics in Acetonitrile solvent molecules surrounding the sites involved in the transfer as well as their orientation in the solvation shell. (ii) The solvation response that we obtain is in excellent agreement with the experimental data13 for the time-dependent Stokes shift within the experimental uncertainties. The model that we use is therefore suitable for an accurate description of time-dependent solute-solvent electrostatic interactions. We have also provided an analysis based on the comparison with a simulation in which we have constrained the C153 geometry to show the improvement that our approach can lead to. According to our result, the short-time dynamics is dominated by the solvent response, while the effect of the solute motion influences the longer time scale. (iii) Even though the potential used in the simulation is not polarizable, the inclusion of the internal degrees of freedom of the solute has allowed us to take into account some effects connected to polarizability. The solute molecular dipole shows a Gaussian-like distribution, which is a similar result to the one obtained in simulations of polarizable, but rigid, C153 in solution.44 We conclude by outlining some of the ways in which our modeling of solvation could be improved. Our model in the present work has not included mutual polarization between the solvent and the excited-state solute. The use of the fluctuating charge model82 in the simulation represents a possible way to extend the approach that we have developed to the inclusion of polarizability effects. A combined investigation between this approach and QM calculations of the ground- and excited-state solute polarizabilities in the presence of the solvent is presently under way. Given that we have shown that intramolecular solute dynamics makes a nonnegligible contribution to S(t), it would be worth including in the model the changes in the solute intramolecular potential associated with the S0 f S1 electronic transition. QM calculations designed to provide us with the appropriate S1 state parameters are planned. References and Notes (1) Simon, J. D. Acc. Chem. Res. 1988, 21, 128. (2) Barbara, P. F.; Jarzeba, W. AdV. Photochem. 1990, 15, 1. (3) Bagchi, B. Annu. ReV. Phys. Chem. 1989, 40, 115. (4) Maroncelli, M. J. Mol. Liq. 1993, 57, 1. (5) Fleming, G. R.; Cho, M. Annu. ReV. Phys. Chem. 1996, 47, 109. (6) Stratt, R. M.; Maroncelli, M. J. Phys. Chem. 1996, 100, 12981. (7) Ladanyi, B. M. Theoretical Methods in Condensed Phase Chemistry. In Theoretical Methods in Condensed Phase Chemistry; Schwartz, S. D., Ed.; Kluwer: Dordrecht, The Netherlands, 2000. (8) Rosenthal, S. J.; Xie, X. L.; Du, M.; Fleming, G. R. J. Chem. Phys. 1991, 95, 4715. (9) Maroncelli, M.; Kumar, P. V.; Papazyan, A.; Horng, M. L.; Rosenthal, S. J.; Fleming, G. R. Ultrafast reaction dynamics and solvent effects. In Ultrafast Reaction Dynamics and SolVent Effects; Gauduel, Y.; Rossky, P. J., Eds.; American Institute of Physics: New York, 1994. (10) Jimenez, R.; Fleming, G. R.; Kumar, P. V.; Maroncelli, M. Nature 1994, 369, 471. (11) Castner, E. W.; Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1987, 86, 1090. (12) Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1987, 86, 6221. (13) Horng, M. L.; Gardecki, J. A.; Papazyan, A.; Maroncelli, M. J. Phys. Chem. 1995, 99, 17311. (14) Gardecki, J.; Horng, M. L.; Papazyan, A.; Maroncelli, M. J. Mol. Liq. 1995, 65-66, 49. (15) Lewis, J. E.; Maroncelli, M. Chem. Phys. Lett. 1998, 282, 197. (16) Gardecki, J. A.; Maroncelli, M. Chem. Phys. Lett. 1998, 301, 571. (17) Luther, B. M.; Kimmel, J. R.; Levinger, N. E. J. Chem. Phys. 2002, 116, 3370. (18) Kovalenko, S. A.; Ruthmann, J.; Ernsting, N. P. Chem. Phys. Lett. 1997, 271, 40. (19) Ruthmann, J.; Kovalenko, S. A.; Ernsting, N. P.; Ouw, D. J. Chem. Phys. 1998, 109, 5466. (20) vanderMeulen, P.; Zhang, H.; Jonkman, A. M.; Glasbeek, M. J. Phys. Chem. 1996, 100, 5367.

J. Phys. Chem. B, Vol. 109, No. 8, 2005 3563 (21) Passino, S. A.; Nagasawa, Y.; Fleming, G. R. J. Chem. Phys. 1997, 107, 6094. (22) Passino, S. A.; Nagasawa, Y.; Joo, T.; Fleming, G. R. J. Phys. Chem. 1997, 101, 725. (23) Lee, S. H.; Lee, J. H.; Joo, T. J. Chem. Phys. 1999, 110, 10969. (24) Mente, S. R.; Maroncelli, M. J. Phys. Chem. B 1999, 103, 7704. (25) Diraison, M.; Millie, P.; Pommeret, S.; Gustavsson, T.; Mialocq, J. C. Chem. Phys. Lett. 1998, 282, 152. (26) Lobaugh, J.; Rossky, P. J. J. Phys. Chem. A 2000, 104, 899. (27) Friedman, H. L.; Raineri, F. O.; Hirata, F.; Perng, B. C. J. Stat. Phys. 1995, 78, 239. (28) Friedman, H. L.; Raineri, F. O.; Perng, B. C. J. Mol. Liq. 1995, 65-66, 7. (29) Ishida, T.; Hirata, F.; Kato, S. J. Chem. Phys. 1999, 110, 11423. (30) Raineri, F. O.; Perng, B. C.; Friedman, H. L. Chem. Phys. 1994, 183, 187. (31) Biswas, R.; Bagchi, B. J. Phys. Chem. A 1999, 103, 2495. (32) Roy, S.; Komath, S.; Bagchi, B. J. Chem. Phys. 1993, 99, 3139. (33) Nishiyama, K.; Yamaguchi, T.; Hirata, F.; Okada, T. Pure Appl. Chem. 2004, 76, 71. (34) Maroncelli, M. J. Chem. Phys. 1991, 94, 2084. (35) Kumar, P. V.; Maroncelli, M. J. Chem. Phys. 1995, 103, 3038. (36) Lobaugh, J.; Rossky, P. J. J. Phys. Chem. A 1999, 103, 9432. (37) Sudholt, W.; Staib, A.; Sobolewski, A. L.; Domcke, W. Phys. Chem. Chem. Phys. 2000, 2, 4341. (38) Ladanyi, B. L.; Perng, B. C. J. Phys. Chem. A 2002, 106, 6922. (39) Ladanyi, B. M.; Stratt, R. M. J. Phys. Chem. 1995, 99, 2502. (40) Ladanyi, B. M.; Stratt, R. M. J. Phys. Chem. 1996, 100, 1266. (41) Ladanyi, B. M.; Klein, S. J. Chem. Phys. 1996, 105, 1552. (42) Ladanyi, B. M.; Maroncelli, M. J. Chem. Phys. 1998, 109, 3204. (43) Cichos, F.; Brown, R.; Bopp, P. A. J. Chem. Phys. 2001, 114, 6824. (44) Cichos, F.; Brown, R.; Bopp, P. A. J. Chem. Phys. 2001, 114, 6834. (45) Reynolds, L.; Gardecki, J. A.; Frankland, S. J. V.; Horng, M. L.; Maroncelli, M. J. Phys. Chem. 1996, 100, 10337. (46) Cichos, F.; Willert, A.; Rempel, U.; VonBorczyskowski, C. J. Phys. Chem. A 1997, 101, 8197. (47) Molotsky, T.; Huppert, D. J. Phys. Chem. A 2002, 106, 8525. (48) Martins, L. R.; Skaf, M. S. Chem. Phys. Lett. 2003, 370, 683. (49) Cance`s, E.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032. (50) Miertusˇ, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (51) Flory, W. C.; Blanchard, G. J. Appl. Spectrosc. 1998, 52, 82. (52) Jones, G. J., II.; Jackson, W. R.; Choi, C.; Bengmark, W. R. J. Chem. Phys. 1989, 89, 294. (53) Mu¨hlpfordt, A.; Schanz, R.; Ernsting, N. P.; Fartzinov, V.; Grimme, S. Phys. Chem. Chem. Phys. 1999, 1, 3209. (54) Cave, R. J.; Castner, E. W., Jr. J. Phys. Chem. A 2002, 106, 12117. (55) Ingrosso, F.; Menucci, B.; Tomasi, J. J. Mol. Liq. 2003, 108, 21. (56) Mennucci, B.; Cammi, R.; Tomasi, J. J. Chem. Phys. 1998, 109, 2798. (57) Carter, E. A.; Hynes, J. T. J. Chem. Phys. 1991, 94, 5961. (58) Hansen, J. P.; McDonald, I. R. Theory of simple liquids; Academic Press: London, 1976. (59) Stell, G.; Patey, G. N.; Høye, J. S. AdV. Chem. Phys. 2981, 48, 183. (60) Edwards, D. M. F.; Madden, P. A.; McDonald, I. R. Mol. Phys. 1983, 51, 1141. (61) Becke, A. J. Chem. Phys. 1993, 98, 5648. (62) Cammi, R.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 110, 9877. (63) Bondi, A. J. Phys. Chem. 1964, 68, 441. (64) Cammi, R.; Cossi, M.; Tomasi, J. J. Phys. Chem. 1996, 100, 4611. (65) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.03; Gaussian, Inc.: Pittsburgh, PA, 2003. (66) Cammi, R.; Mennucci, B.; Tomasi, J. J. Phys. Chem. A 2000, 104, 5631. (67) Smith, W.; Forester, T. R. DL_POLY_2; CCP5: Daresbury, Warrington, England, 1999.

3564 J. Phys. Chem. B, Vol. 109, No. 8, 2005 (68) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1989. (69) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (70) Gallant, R. W. Hydrocarbon Process. 1969, 48, 135. (71) Bo¨hm, H. J.; McDonald, I. R.; Madden, P. A. Mol. Phys. 1983, 49, 437. (72) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (73) Watkins, E. K.; Jorgensen, W. L. J. Phys. Chem. A 2001, 105, 4118. (74) Rizzo, R. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1999, 121, 4827.

Ingrosso et al. (75) Nugent, S.; Ladanyi, B. M. J. Chem. Phys. 2004, 120, 1. (76) Egorov, S. A. J. Chem. Phys. 2000, 113, 1950. (77) Adams, J. E. J. Phys. Chem. B 1998, 102, 7455. (78) Biswas, R.; Lewis, J. E.; Maroncelli, M. Chem. Phys. Lett. 1999, 310, 485. (79) Ando, K.; Kato, S. J. Chem. Phys. 1991, 95, 5966. (80) Fonseca, T.; Ladanyi, B. M. J. Phys. Chem. 1991, 95, 2116. (81) Aherne, D.; Tran, V.; Schwartz, B. J. J. Phys. Chem. B 2000, 104, 5382. (82) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141.