J. Phys. Chem. A 2004, 108, 7347-7355
7347
Simulations of the Vibrational Relaxation of a Model Diatomic Molecule in a Nanoconfined Polar Solvent Shenmin Li,† Tricia D. Shepherd,‡ and Ward H. Thompson*,† Department of Chemistry, UniVersity of Kansas, Lawrence, Kansas 66045-7582, and Department of Chemistry/ Physics, Westminster College, Salt Lake City, Utah 84105 ReceiVed: April 14, 2004; In Final Form: June 25, 2004
The vibrational dynamics of a model diatomic anion solute dissolved in a methyl iodide solvent confined in a nanoscale spherical cavity are investigated by molecular dynamics simulations. The effect of confining the solvent on the vibrational energy relaxation time T1, solvent-induced frequency shift 〈δω〉, and pure dephasing time T2* is examined by comparing the results from confined systems of varying size (cavity radius 0.8-2 nm) to those from the bulk system. It is found that T1 increases monotonically toward the bulk solvent value with increasing cavity size, and good agreement is found between equilibrium molecular dynamics simulations based on perturbation theory and classical nonequilibrium dynamics simulations. In contrast to T1, the solventinduced frequency shift and the dephasing time do not change monotonically with cavity size. The results are discussed in terms of the changes in solvent structure and dynamics due to confinement.
1. Introduction The structure and dynamics of confined liquids have attracted increasing attention in recent years.1-14 The confinement of liquids in microporous media, such as sol-gels, zeolites, supramolecular assemblies, and reverse micelles, can lead to dramatic changes in their static and dynamical properties from those of the bulk system.15 A number of experimental approaches have been used to probe these changes including steady-state absorption and fluorescence spectra,16-33 timedependent fluorescence,16-33 NMR spectroscopy,1,34-37 OKE spectroscopy,38-43 dielectric relaxation spectroscopy,44,45 chargetransfer reactions,46,47 Raman spectroscopy,2,3,9,10,48 and infrared pump-probe experiments.6,7 Taken together these experiments can access equilibrium and nonequilibrium structure and dynamics including processes dominated by long- and shortrange interactions. Among these, the vibrational dynamics investigated by Raman spectra and time-resolved vibrational pump-probe experiments have been comparatively less wellstudied. The results of these measurements are vibrational frequency shifts (〈δω〉), dephasing times (T2*), and energy relaxation times (T1), processes in which the solvent-solute short-range repulsive interactions play a large role. Thus, they provide a valuable complement to measurements that are sensitive to the longer-range electrostatic interactions (e.g., optical spectra); these quantities also each probe slightly different solvent properties and solute-solvent interactions. In addition, they are a useful starting point for well-controlled studies of energy transfer across nanoscale interfaces, a topic of increasing interest. In this paper we investigate the vibrational dynamics in a nanoconfined solvent by calculating 〈δω〉, T2*, and T1 by molecular dynamics simulations. Recently, Zhong et al. have studied vibrational energy relaxation (VER) of pseudohalide ions, N3-, NCO-, and NCSconfined in nanoscale water pools inside nonionic reverse micelles (RMs) by ultrafast infrared spectroscopy.6,7 These three † ‡
University of Kansas. Westminster College.
systems all gave similar results as a function of the RM water pool radius; i.e., the VER times are about three times longer for the smallest RM studied than the values measured for bulk water, and become shorter with increasing RM size (but do not reach the bulk value for the largest RM studied). The solventinduced frequency redshift was also observed to decrease toward the value in bulk water as the RM size was increased.5 However, the magnitudes of the frequency shifts are much different for various solutes or even the same solute in different micelles. For instance, the redshifts for NCS- are much smaller than those for NCO- and N3-. The redshift is less than one wavenumber for NCS- in bulk H2O but in nonionic RMs it is 7-8 cm-1. Zhong et al. qualitatively explained the longer VER times and solvent-induced redshifts in terms of “the reduced interactions between the ions and solvent” which is ultimately related to the water-surfactant interactions.7 Interestingly, several groups1-3,9,10,48,49 have measured the vibrational dephasing time, T2, of neat, nonaqueous solvents in sol-gels and found a different trend. Specifically, Jonas and co-workers found that T2 becomes shorter with decreasing pore size (in fact, T2-1 ∝ 1/Rcav), and the solvent induces a frequency blueshift. The accelerated vibrational dephasing upon confinement was attributed to an increase in the orientational order relative to the bulk liquid and to liquid-surface interactions. Tominaga et al. obtained similar results.49 Both groups explained their observations in terms of a two-state model in which liquid molecules inside the pore are classified into “surface molecules” and “bulk molecules”.1,48,49 These experimental results indicate some of the possible interesting vibrational dynamics occurring in confined solvents and demonstrate the sensitivity to the characteristics of the confining framework. However, there has been little theoretical study of confinement effects on vibrational dynamics.4,50-52 In this paper, we investigate the vibrational dynamics in a simple model system of a solute dissolved in a nanoconfined solvent. The emphasis here is on developing a general understanding of the origins of the trends in 〈δω〉, T2*, and T1 for a simple model system. The present work is not intended as a realistic simulation
10.1021/jp048361e CCC: $27.50 © 2004 American Chemical Society Published on Web 08/13/2004
7348 J. Phys. Chem. A, Vol. 108, No. 36, 2004
Li et al.
of a system studied experimentally, which will be a focus of future work. Thus, the solute is taken to be a model diatomic anion (AB-), and the solvent is dipolar, non-hydrogen bonding CH3I. In this initial study we choose a confinement framework that is relatively inert: a rigid, hydrophobic spherical cavity that can have, at best, a limited participation in the vibrational dynamics. This serves two purposes: (1) it allows for comparisons with previous simulations of steady-state absorption and fluorescence spectra4 and time-dependent fluorescence measurements, and (2) it provides a useful reference point for investigating the role of the surface. The vibrational energy relaxation time, solvent-induced frequency shift, and pure dephasing time are calculated for spherical cavities of varying size by equilibrium molecular dynamics simulations (EMD) and compared to the values in the bulk system. As an important validation step, nonequilibrium molecular dynamics simulations (NEMD) are also used to calculate the VER time. The rest of this paper is organized as follows. In section II, we describe the basic theory of vibrational energy relaxation and vibrational dephasing. A brief discussion of the method used for computing the numerical Fourier transform is also given in this section. The details of the computational model and procedures are given in section III. The calculated results are presented and possible vibrational relaxation mechanisms in the confined system are discussed in section IV. Finally, the main conclusions are summarized in section V. II. Theoretical Background A. Vibrational Energy Relaxation. 1. Equilibrium Molecular Dynamics (EMD). The most prevalent approach for calculating the vibrational energy relaxation time T1 is based on perturbation theory53-55
T1 )
µkBT ξ(ω0)
(1)
where µ is the reduced mass of the solute molecule, kB is Boltzmann’s constant, T is the temperature, and ξ(ω0) is the power spectrum of the force-force time correlation function, 〈δF1(t)δF1(0)〉, at the solute vibrational frequency ω0:
ξ(ω0) )
∫0∞ cos(ω0t)〈δF1(t)δF1(0)〉 dt
(2)
Here
[
F1 ) µ
]
FA B FB B L2 ‚rˆAB + mA mB Ire
(3)
where B FA and B FB denote the total forces on the A and B atoms, respectively, rˆAB is the unit vector from A to B, L is the solute angular momentum, I is the of moment of inertia, and re is the equilibrium AB bond distance. Thus, δF1(t) is the fluctuation of the force F1 exerted by solvent atoms and solute rotation along the solute bond at time t. In principle, the force TCF should be evaluated quantum mechanically. However, since this is computationally unfeasible for a many-body real system, the classical force TCF is usually used and a correction factor is introduced to account for quantum effects.56-59 The classical TCF can be easily obtained from the EMD simulations, by monitoring the fluctuations of force F1 at each MD time step. There are well-known issues concerning the accuracy of the EMD approach for calculating T1. The first is how to deal with the Fourier transform of the force-force time correlation function (TCF), particularly for a solute with a high vibrational frequency.
The second is how to properly account for quantum effects, which has focused attention on the calculation of quantum correction factors (QCFs) to the classical force TCF. Two kinds of techniques are widely used to perform the Fourier transform. One is fitting the force TCF to an analytic function from which an analytic Fourier transform form can be obtained.60,61 This method is usually used for short-time fitting and therefore is primarily applicable when the force TCF decays quickly. The other is performing the Fourier transform numerically. However, the traditional fast Fourier transform can suffer from numerical noise for solutes with high vibrational frequency.62,63 In this paper, we use the modified Hurwitz-Zweifel (MHZ) method64,65 to calculate the cosine Fourier integral in eq 2 numerically. This method was adopted to calculate the VER time of I2 in a xenon solvent and yielded satisfactory results.62 With regard to the QCFs,56-59 we do not use them in our calculations for two reasons: (1) there is not an obvious form of the QCF to use for an arbitrary liquid system,59 and (2) our interest is in the confinement effect on T1 and there is not yet a clear method for estimating the dependence of the quantum correction to the VER time on the confined system size. The MHZ approach, presented by Thakkar and Smith,65 is to subdivide the range and integrate between the successive zeros of cos(ω0t), thus converting the infinite integral to a summation of terms Cn, which is referred to as the Hurwitz-Zweifel expansion.64 The Clenshaw-Curtis quadrature method66 is used here to evaluate the Cn as suggested by Thakkar and Smith.65 The main difficulty is that the Cn series may converge slowly. To accelerate the convergence, a reliable and fairly general scheme is to use the van Wijngaarden modification of the Euler transformation.67 The details of the numerical evaluation procedures can be found in reference 65. 2. Nonequilibrium Molecular Dynamics (NEMD). Nonequilibrium molecular dynamics can be used to determine classical vibrational energy relaxation times directly by investigating the time decay of excess vibrational energy. If the energy decay is exponential, the relaxation time T1 is given by
〈EV(t)〉 - 〈EV(∞)〉 〈EV(0)〉 - 〈EV(∞)〉
) e-t/T1
(4)
where 〈Eν(t)〉 is the nonequilibrium average of the vibrational energy for a diatomic solute at time t, and 〈Eν(∞)〉 and 〈Eν(0)〉 are the equilibrium average vibrational energy and the nonequilibrium average initial energy, respectively. The solute vibrational energy at time t is monitored in a nonequilibrium trajectory such as
EV(t) )
pr2(t) + V[r(t)] 2µ
(5)
where r is the solute anion bond distance with conjugate momentum pr, and V(r) is the solute potential energy. B. Vibrational Dephasing and Frequency Shifts. Theoretical treatments have show that the dephasing time T2 and VER time T1 are related by68
1 1 1 ) + T2 2T1 T2*
(6)
where T2* is the pure dephasing time. It is generally the case in a liquid that T1 is comparatively long and makes a negligible contribution to the overall vibrational dephasing time T2. Therefore, in many cases the spectral line width is a measure of the pure dephasing time T2* only.2,3 According to the Kubo
Simulations of Vibrational Relaxation
J. Phys. Chem. A, Vol. 108, No. 36, 2004 7349
theory69-71 of dephasing, the pure dephasing time T2* is given by
1 ) T2*
∫0
∞
〈δω(0)δω(t)〉 dt
(7)
where 〈δω(0)δω(t)〉 is the time correlation function of solventinduced vibrational frequency shift, δω(t) ) ω(t) - 〈ω〉. A widely used method for calculating δω(t) in EMD simulations was introduced by Oxtoby.72 Using time-independent perturbation theory and considering only the first excited state, the frequency shift can be described by a linear term (δω1) and a quadratic term (δω2)
δω(t) ) δω1(t) + δω2(t)
(8)
TABLE 1: Potential Parametersa site
(kcal/mol)
(9)
δω2(t) ) [(Q2)11 - (Q2)00]F2(t)
(10)
Here, F1(t) is the force defined in eq 3 while F2(t) is half the derivative of F1(t) with respect to the solute vibrational coordinate displacement, Q ) r - re, evaluated at Q ) 0. In this work, F2(t) is calculated numerically during the equilibrium MD simulation. The expectation values Q11, Q00, (Q2)11, and (Q2)00 are calculated by solving the stationary Schro¨dinger equation of the isolated solute molecule, H ˆ ψ ) Eψ, where the Hamiltonian has the form
H ˆ )
pˆ µ2 + Vˆ (Q) 2µ
(11)
and the diatomic potential V is taken to be a Morse function. In addition, for a rigid solute diatomic molecule, the vibrationrotation coupling contribution to the frequency shift can be split as in eq 873
δωrot(t) ) δωr1(t) + δωr2(t)
(12)
A B
0.1195 0.1195
CH3 I
0.2378 0.5985
3.77 3.83
0.46
Cavity Wall 2.5 0.0
a
δωr2(t) ) 3[(Q )11 - (Q )00] 2
2
(13)
(14)
2Ire2
where L(t) is the solute angular momentum at time t, and since 〈L2〉/2I ) kT, this gives
{
}
(Q11 - Q00) (Q2)11 - (Q2)00 +3 kT (15) 〈δωrot〉 ) -2 re r2 e
This implies that within perturbation theory the frequency shift induced by solute vibration-rotation coupling is independent of the surroundings. III. Model and Simulation Details Molecular dynamics simulations are carried out to investigate the vibrational dynamics of one solute AB- anion in a CH3I solvent confined in spherical cavities of varying size. The corresponding results for the bulk system are also presented. All simulations are carried out at an average temperature of
1.2
15 15
2.16
15 126.8
T1-EMD† (ps)a T1-EMD (ps) T1-NEMD (ps) T2* (ps) no. of solvent molecules
Rcav ) 8Å
Rcav ) 10 Å
Rcav ) 12 Å
Rcav ) 15 Å
Rcav ) 20 Å
bulk
47.0 33.5 37.0 0.80 11
54.0 44.3 39.5 0.04 24
58.7 46.9 43.0 0.13 44
61.1 51.9 48.6 0.51 92
62.9 53.0 49.1 0.71 234
82.3 63.5 63.5 3.67 107
a T -EMD† denotes the VER time obtained by EMD simulation 1 without the solute vibration-rotation coupling contribution.
298 K and a solvent density of 2.0 g/cm3. For simulations of confined solvents, the cavity radius is chosen to be 8, 10, 12, 15, and 20 Å; the corresponding numbers of solvent molecules at this density are given in Table 2. For simulations of the bulk system, the minimum image convention and periodic boundary conditions are used, with 107 CH3I molecules in a cubic box of length 23.29 Å. The form of the interaction potential is the same as in our previous work.4 The potential parameters are given in Table 1. For the confined system, the potential functions are a sum of Coulomb, uC(rij), and Lennard-Jones, uLJ(rij), interactions for all pairs of sites and molecule-wall interaction terms, uw(ri), for all sites:
uC(rij) + ∑uLJ(rij) + ∑uw(ri) ∑ i
