Molecular Dynamics Study of the Temperature-Dependent Optical

J. Phys. Chem. B 2008, 112, 7837–7849

7837

Molecular Dynamics Study of the Temperature-Dependent Optical Kerr Effect Spectra and Intermolecular Dynamics of Room Temperature Ionic Liquid 1-Methoxyethylpyridinium Dicyanoamide Zhonghan Hu,† Xuhui Huang,‡ Harsha V. R. Annapureddy,† and Claudio J. Margulis*,† Department of Chemistry, The UniVersity of Iowa, Iowa City, Iowa 52242, and Department of Chemistry, Columbia UniVersity, New York, New York 10027 ReceiVed: January 24, 2008; ReVised Manuscript ReceiVed: March 11, 2008

We have performed classical molecular dynamics simulations to calculate the Optical Kerr effect (OKE) spectra of 1-methoxyethylpyridinium dicyanoamide, a room-temperature ionic liquid (IL) which has been recently studied by Shirota and Castner (Shirota, H.; Castner, E. J. Phys. Chem. A 2005, 109, 9388-9392) in comparison to its neutral isoelectronic solvent mixture. Our theoretical and computational studies show that the decay of the collective polarizability anisotropy correlation exhibits several different time scales originating from inter- and intramolecular dynamics, in good agreement with experiments. What’s more, we find that the portion of the collective anisotropic polarizability relaxation due to“interaction-induced” phenomena is important at times much longer than those observed in normal solvents when these are far from their glass transition temperature. From our long (60 ns) molecular dynamics simulations, we are able to determine the appropriate time scales for orientational relaxation and interaction-induced processes occurring in the liquid. We find that the cationic contribution to the OKE signal is predominant. Because of the slow nature of relaxation processes in ILs, these calculations are very time, memory, and storage intensive. In the context of this research, we have developed a polarizable force field for this system and also theoretical methodology to generate molecular polarizabilities for arbitrarily shaped molecules and ions from corresponding atomic polarizabilities. We expect this methodology to have an important impact on the speed of molecular dynamics simulations of polarizable systems in the future. 1. Introduction The intermolecular dynamics of room-temperature ionic liquids (RTILs) has been the subject of many recent theoretical and experimental studies; see, for example, refs 1–28 and citations therein. This list is by no means comprehensive. In this article, we focus our attention on one particular technique, the Optical Kerr Effect spectroscopy, which can be used to probe intramolecular as well as collective intermolecular dynamics of liquids.6,8,9,11,18,24,29,30 For complex solvents, the experimental analysis of the lowfrequency part of the Kerr spectrum is problematic and calls for the use of molecular dynamics simulations to elucidate the different dynamical processes that give rise to the observed signal. At very low frequencies, the spectra are usually dominated by rotational diffusive processes, while slightly higher frequencies usually correspond to librational and, in general, intermolecular vibrational motions. The spectra at higher frequencies report of intramolecular vibrational modes that are Raman-active in the liquid. The OKE spectrum can be derived from a polarizabilitypolarizability correlation function. At long times,8 to a first approximation, this correlation function may be considered to be similar to a reorientational correlation function. Keeping the terminology from scattering theory,9 at shorter times, inertial and librational motion are usually termed “collision-induced dynamics”. As Hunt and co-workers describe in their interesting article,31 fast subpicosecond processes are better studied in the * To whom correspondence should be addressed. † The University of Iowa. ‡ Columbia University.

frequency domain, while longer time response processes can be better analyzed in the time domain. If the dynamics is followed for long enough times (i.e., long pulse delays), the temperature dependence of the OKE spectrum can report on the nature of the rotational relaxation in ionic liquids. For example, Cang and co-workers8 have been able to fit their timeand temperature-dependent data to a sum of power laws multiplied by an exponential decay. From the time constant in the exponential, they have argued that their systems satisfy the Debye-Stokes-Einstein (DSE) model for rotational diffusion. It is well understood that an increase in the strength of intermolecular forces in liquids is strongly correlated with an increase in the complexity of the low-frequency part of the OKE spectrum.9 As an example, librations in ionic liquids appear at higher frequencies that those in DMSO, and a combination of several functions is required to fit the low-frequency part of the OKE spectrum in the case of ILs, which reflects the inhomogeneous nature of these systems. In systems with strong Coulomb interactions such as in ILs, the molecular polarizability of the ions should be strongly influenced by the dynamics of neighboring molecules. This is in contrast to the case in many other solvents in which the change in polarizability after initial inertial motion is solely dominated by the trivial change in polarizability observed in a fixed laboratory frame due to molecular rotations. The degree to which rotations and other collision-induced or interaction-induced processes contribute to the polarization relaxation is difficult to derive from an analysis of the spectroscopic data, but computer simulations and the use of projection operators can help enormously in separating

10.1021/jp800729g CCC: $40.75  2008 American Chemical Society Published on Web 06/07/2008

7838 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Figure 1. Cation [MOEPY+] (A) and anion [DCA-] (B).

contributions and defining time scales on which different processes are active. Perhaps one of the most important problems commonly found in the literature is that OKE results are only reported on the subnanosecond or several nanosecond time scales. Computer simulations have shown that for many IL systems, diffusive rotational and translational dynamics is only established at room temperature on a longer time scale. An interesting phenomenon reported by Quitevis and co-workers6,11,24 is that mixtures of ILs appear to show OKE spectra that are additive. These findings are reminiscent of other studies by the group of Professor Hamaguchi32,33 in which Raman signals were found to be reminiscent of those of distinct solid polymorphs in the liquid phase. Shirota and Castner34 have carried out interesting OKE studies on low-viscosity silicon-based ILs. Particularly interesting to us is their study on an ionic liquid and its isoelectronic molecular mixture counterpart. We have chosen to focus our attention on the system studied by Shirota and Castner, 1-methyoxyethylpyridinium dicyanoamide ([MOEPY+][DCA-]; see Figure 1), because their detailed time and frequency domain spectra provided us with a great starting point for attempting a comprehensive understanding of the intermolecular vibrational and diffusive dynamics in ILs. Shirota and Castner demonstrated that the low-frequency Kerr spectrum of the IL appears at a higher frequency than its neutral counterpart and also that the signal has a much broader width. In order to describe the decay of the Kerr signal, Shirota and Castner chose multiple exponential functions. The slowest time constant obtained from the analysis of their data in the case of the IL was 1.08 ns, while in the case of the neutral mixture, this constant was 29 ps. It is very important to notice however that a full relaxation of the Kerr intensity could not be observed within the 750 ps time window of their experiment, leaving open the possibility that processes with even longer time scales could be present in the liquid. As far as we know, no theoretical studies have directly addressed the time scales for intermolecular interactions and collective reorientational dynamics by strict comparison and calculation of the OKE transient intensities in the area of roomtemperature ionic liquids. In part, this is because of the tremendous computational difficulties that must be overcome to produce reliable long time data in the case of RTILs. First, there is the issue of long time rotational relaxation, which requires extremely long simulations (on the order of hundreds of nanoseconds for [MOEPY+][DCA-] at room temperature, which corresponds to many months of parallel simulations on powerful computer clusters) in order to fully converge the polarization correlation function required to compute the OKE transient signal. A second difficulty arises from the lack of accurate polarization parameters for ionic liquids (in fact, the whole area of force field development for ILs is in its infancy when compared, for example, to that for water models, though some good polarizable force fields have been developed in the Voth group62,63). The third difficulty is the time, RAM, and disk storage required for the calculation of collective polarizabilities along simulations based on the well-known dipole-induceddipole model.35–38 Other than using faster computers and hoping for more experimental structural data in the form of neutron scattering, there is not much that we can do to alleviate issues

Hu et al. one and two. The third issue, on the contrary, lends itself to the development of theoretical tools that can reduce the cost of working with atomic polarizabilities. In this context, we have developed an algorithm of general applicability to generate approximate molecular polarizability parameters for arbitrarily shaped molecules or ions from corresponding atomic polarizabilities. This method tremendously reduces memory requirements, allowing for fast inversion of much smaller symmetric matrices and can therefore be used to study larger systems with very little decrease in accuracy when compared to the widely used self-consistent36,38 and inverse matrix algorithm39,40 to solve the DID model. The atomic polarizabilities that we transform using our algorithm into molecular polarizabilities are those generated by the methodology developed by Stern et. al.41 Due to the arduous challenges described above, our simulations were performed at higher temperatures (400, 500, and 600K) than those studied by Shirota and Castner, and our system is somewhat small (64 pairs of ions). Despite this, we are able to capture virtually all of the features of the OKE frequency domain and time domain transient spectra observed experimentally. From our calculations, we are able to obtain the time scales and the degrees to which different rotational and collisioninduced components contribute to the signal decay as well as the anionic and cationic contributions to the overall signal. This paper is organized as follows. For completeness, in section 2, we will present the basic equations involved in OKE spectroscopy and will describe the standard DID model used in the past for the calculation of the OKE spectra in traditional solvents by other research groups.36,38,40 In section 3, we will describe our algorithm to calculate molecular polarizabilities from atomic parameters and will show a comparison between our results and those obtained using as a gold standard the selfconsistent method with full atomic polarizabilities. Molecular dynamics simulation details for the [MOEPY+][DCA-] system are provided in section 4. Atomic and molecular polarizabilities as well as other force field parameters are provided in the Appendix. Finally, we analyze our results in light of the experiments by Shirota and Castner in section 5 and draw conclusions in section 6. 2. Theoretical Background 2.1. Optical Kerr Effect. The Optical Kerr Effect experiment measures the transient anisotropy in the refractive index ∆n(t)36,42

∆n(t) )

∫-∞t dτ Ipump(τ)R(t - τ)

where R(t) is the third-order nonlinear response function, which has an instantaneous hyperpolarizability part (arising from the electronic response) and a time correlation part

R(t) ) ζδ(t) - C˙2(t) C˙2(t) is the first derivative of the anisotropic collective polarizability-time correlation function. The observed heterodynedetected transmission of a Kerr cell is approximately given by

IOKE(t) )

∫-∞t dτ G0(τ)R(t - τ)

where G0(t) is the zero background laser pulse intensity autocorrelation

G0(t) )

∫-∞∞ dτ Ipump(t)Iprobe(t - τ)

Because the third-order response function R(t) is a combination of a hyperpolarizability ζδ(t) and a time correlation function

MD Study of 1-Methoxyethylpyridinium Dicyanoamide

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7839

C˙2(t), the instantaneous response in the OKE transient signal resembles the blank convolution signal G0(t). For clarity, experimentalists almost always subtract this instantaneous contribution when presenting their transient intensities. Our simulations only track the transient signal from the retarded time ˙ 2(t), which probes the inter- and intramocorrelation function C lecular dynamics of the liquid. We directly write this part of the OKE signal as

IOKE(t) ∝ -

∫-∞ dτ G0(τ)C˙2(t) t

j) D j - 1 Tr(D j )jI K 3

(

The anisotropic component of the total collective polarizability j is the sum of that of the molecular polarizablity and that of Π the interaction-induced polarizability

j )Q j +K j Π 2 j can be written as Thus, the time correlation function of Π

j (0) · Π j (t)] 〉 C2(t) ) 〈Tr[Π 2 2 j j j j (t)] 〉 + 〈 〉 〈 [ ] [ ) Tr Q(0) · Q(t) + Tr K(0) · K

(1)

In our calculations, the shape of the probe and pump laser pulses is taken to be a Gaussian function with a 37 fs full width at ˙ 2(t), half-maximum.18 Denoting the Fourier transform of G0(t), C and IOKE as G0(ω), C˙2(ω), and IOKE(ω), we have43

(2)

The above deconvoluted spectrum C˙2(ω) is usually much smoother than the one obtained from a direct Fourier transform of C˙2(t). 2.2. TCF of the Collective Polarizability. The collective j of a group of molecules is the sum of polarizability Π polarizabilities for all isolated molecules. Each isolated molecular polarizablity seen from a laboratory-fixed set of Cartesian coordinates is modified by self-rotation and by interactions of the molecule (or molecular ion) with its surroundings. The standard point dipole/induced point dipole (DID) model expresses the isolated polarizability π j (j) for molecule j as a sum of a molecular component (gas-phase polarizability) and an interaction-induced component N

∑ Tj jk · πj (k)

(3)

k*j

where the dipole-dipole tensor between molecules j and k can be written as

j ) 1 (3rˆrˆ - jI )| T jk r)rjk r3 and rˆ ) r/r is the unit vector along the direction r. Equation 3 can be solved self-consistently36 by taking π j (k) ) R j (k) as the j for a system of N ansatz. The total collective polarizability Π particles is

j ) Π jI + Π j ) Π 0 2

N

∑ πj (j)

)

N

N

j

j

k*j

∑ Rj(j) + ∑ ∑ Rj(j) · Tj jk · πj (k)

j +D j )M 1 j )jI + M j - 1 Tr(M j )jI + 1 Tr(D j )jI + D j - 1 Tr(D j )jI ) Tr(M 3 3 3 3 1 j )jI + Q j + 1 Tr(D j )jI + K j ) Tr(M 3 3

(

)

(

)

where the anisotropic component of the molecular polarizability j and interaction-induced polarizability D j M

j) M j - 1 Tr(M j )jI Q 3

(

and

)

(4)

I MI ) CM 2 (t) + C2(t) + C2 (t)

j ]〉 〈 [j CM 2 (t) ) Tr Q(0) · Q(t) j (0) · K j (t)] 〉 CI2(t) ) 〈Tr[K j j j ]〉 〈 [j CMI 2 (t) ) Tr Q(0) · K(t) + K(0) · Q(t) Equation 4 separates the TCF C2(t) for the anisotropic collective j 2 into three components: a molecular component polarizability Π CM , an interaction-induced component CI2(t), and a molecular 2 interaction-induced cross-correlation component C2MI(t). This separation enables one to distinguish the pure molecular contribution from the pure interaction-induced contribution. The problem with this partitioning is that it really does not say much about the dynamics in the liquid. All components follow, to some extent, the librational and, in general, reorientational dynamics of the molecular ions. In order to separate out the fraction of the interaction-induced component that tracks the dynamics of molecular reorientation, it is informative to project the interactionj onto the molecular anisotropic induced anisotropic polarizability K j ,35,44 which, at low enough frequencies when polarizability Q intramolecular processes are not resonant, is presumed to follow the orientation of the molecule with respect to fixed laboratory axis. The projected component is

j j j ) Tr[K · Q] Q j ∆Q j ·Q j] Tr[Q where the coefficient ∆ is defined as

∆≡

j ·Q j] Tr[K j ·Q j] Tr[Q

j 2 can be Therefore, the anisotropic collective polarizability Π separated into a pure reorientational part and a pure collisioninduced part (without reorientation)

j

N

j (0) · K j (t) + K j (0) · Q j (t)] 〉 〈Tr[Q where

IOKE(ω) C˙2(ω) ) G0(ω)

j (j) + R j (j) · π j (j) ) R

)

jR+Π j CI j )Q j +K j ) (1 + ∆)Q j + (K j - ∆Q j))Π Π 2 2 2 The corresponding time correlation functions are RI C2(t) ) CR2 (t) + CCI 2 (t) + C2 (t)

(5)

j (0) · (1 + ∆(t))Q j (t)] 〉 CR2 (t) ) 〈Tr[(1 + ∆(0))Q j ) (j j )] 〉 〈 [( j CCI 2 (t) ) Tr K(0) - ∆(0)Q(0) · K(t) - ∆(t)Q(t) jR j CI j CI jR CRI 2 (t) ) 〈Tr[Π2 (0) · Π2 (t) + Π2 (0) · Π2 (t)] 〉 RI where CR2 (t), CCI 2 (t), and C2 (t) stand for the reorientational contribution, the collisional-induced contribution, and the cross term contribution to C2(t), respectively.


Hu et al.

The system of interest [MOEPY+][DCA-] can be regarded as a binary mixture of cations and anions. Therefore, we can further decompose the anisotropic component of the total collective polarizability as

j )Π j++Π jΠ 2 2 2 The cationic and anionic components specified in the above equation are

j+) Π 2

∑ πj 2(j+) j+

j-) Π 2

∑ πj 2(j-) j-

where the index j+ (j-) goes over all cationionic (anionic) atoms. Correspondingly, C2(t) can be split into a catonic component, an anionic component, and their cross contribution +/C2(t) ) C+ 2 (t) + C2 (t) + C2 (t)

j+ j+ C+ 2 (t) ) 〈Tr[Π2 (0) · Π2 (t)] 〉 jjC2 (t) ) 〈Tr[Π2 (0) · Π2 (t)] 〉 j+ jjj+ C+/2 (t) ) 〈Tr[Π2 (0) · Π2 (t) + Π2 (0) · Π2 (t)] 〉

(6)

The first derivative of C2(t) can be written as

˙˙ +/C˙2(t) ) C˙+ 2 (t) + C2 (t) + C2 (t)

(7)

Since this function is mostly negative, for convenience, we define D ≡ -C˙2, and we plot this quantity instead.

Equation 3 can be solved self-consistently to obtain the collective polarizability for a group of N atoms under the DID approximation. An alternative procedure is to recast this equation in the form of a 3N by 3N interaction matrix N

∑ Tj jk · πj (k) ) jI k*j

[

matrix.39

and then invert this This equation can be rewritten as follows

j -1(1) R

j -T 12

···

j -T 1N

j -T 21

j -1(2) R

···

j -T 2N

···

···

· · ·

···

j -T N1

j -T N2

···

j -1(N) R

jI

π j (1)I

·

π j (2)I

)

jI

· · I·

···

π j (N)I

jI

where each element of the above 3 × N by 3 × N matrix is a 3 by 3 matrix. The inverse of the 3N by 3N matrix is

[

j j j ··· B B 1N 11 B12 j j j · · · B B B 2N 22 B˜ ) 21 ··· ··· ··· ··· j j j B N1 BN2 · · · BNN

]

j ij of the 3N by 3N matrix B˜ is a 3 by where each component B 3 tensor. The collective polarizability can be written as

N

∑

π j (j) )

j

N

N

i

j

∑ ∑ Bj ij

The matrix B˜ is symmetrical because the gas-phase polarizability j ij are both symmetric. Performing R and dipole-dipole tensor T j can be done the inversion of a symmetrical matrix to get Π quite efficiently, but the memory requirements become very large as the number of atoms increases. In order to facilitate the calculation and reduce the use of memory, we developed the AMPT approximation. This approximation, which dramatically reduces the size of the system of equations by casting the problem in terms of molecular polarizabilities, can be implemented with little loss of accuracy, regardless of whether one decides to solve a self-consistent problem, invert a matrix, or use an extended Lagrangian approach. We have chosen to do a direct matrix diagonalization since, under the AMPT approximation, this was very computationally efficient. The number of atoms in our 64 ion pair simulations is N ) 64 × (22 + 5). This leads to an atomic interaction matrix with 3 × 64 × (22 + 5) ) 5184 by 5184 elements. This interaction matrix can be decomposed into 64 × 64 + 64 × 64 + 2 × 64 × 64 ) 16384 blocks, among which 64 × 64 blocks (each block is a 3 × 22 by 3 × 22 matrix) correspond to cation-cation interactions, another 64 × 64 blocks (each block is a 3 × 5 by 3 × 5 matrix) correspond to anion-anion interactions, and the remaining 2 × 64 × 64 blocks (each block is a 3 × 22 by 3 × 5 matrix) correspond to cation-anion interactions. Under the AMPT approximation, each of these 16384 blocks (irrespective of their size) is replaced by a 3 × 3 molecular interaction matrix. This results in a size reduction from an original 5184 × 5184 atomic interaction matrix into a 3 × 128 by 3 × 128 molecular interaction matrix. The total collective polarizability is obtained by inversion of this molecular interaction matrix. The original inversion method for the DID model expresses the collective polarizability for two interacting molecules a and b as follows

[ ][ ] [ ]

3. The Atomic to Molecular Polarizability Transformation (AMPT) Algorithm

j -1(j) · π R j (j) -

j) Π

j j T a1a2 · · · Ta1aL j j a2a1 Ta2a2 · · · Ta2aL ··· ··· ··· ··· j jT j aLa1 TaLa2 · · · TaLaL -----------j j j T b1a1 Tb1a2 · · · Tb1aL j j j ··· T T T j T a1a1 jT

b2a1

b2a2

b2aL

] [ ][]

j j T a1b2 · · · Ta1bM j j a2b1 Ta2b2 · · · Ta2bM ··· ··· ··· ··· j jT j aLb1 TaLb2 · · · TaLbM -- --------- · j jT j b1b1 Tb1b2 · · · Tb1b2 j j j ··· T T T j T a1b1 jT

b2b1

b2b2

b2b2

··· ··· ··· ··· ··· · ·· ··· ··· j j jT j jT j bMa1 bMa2 · · · TbMaL TbMb1 TbMb2 · · · TbMbM j a1) π( j a2) π(

jI jI

··· j aL) π(

--j b1) π(

··· jI ) -jI

j b2) π( ··· j bM) π(

jI ··· jI

where L and M are the number of atoms for molecules a and b, respectively. From now on, we will refer to the first matrix on the left-hand side as the T interaction matrix. For



Figure 2. Comparison between the accurate atomistic self-consistent approach and molecular-based polarizabilities obtained from our AMPT approach. The time correlation function C2(t) (A) and its negative derivative D(t) (B) are computed from a simulation at T ) 600 K.

j a1a1 ) R convenience, we have used T j -1(1) as the polarizability inverse of the first atom of molecule a and Tj a1b1 ) -Tj 1(L+1) as the interaction between atom 1 from molecule a and atom L + 1 from molecule b. The rest of the matrix element labels have similar meaning. We denote the inverse of the above 3 × (L + M) by 3 × (L + M) T interaction ˜ matrix as H

[

j H a1a1 j H a2a1

j H a1a2 · · · j Ha2a2 · · ·

j H a1aL j H a2aL

j H a1b1 j H a2b1

j H a1b2 · · · j Ha2b2 · · ·

j H a1bM j H a2bM

··· ··· ··· ··· ··· ··· ··· ··· j j j j j j HaLa1 HaLa2 · · · HaLaL HaLb1 H aLb2 · · · HaLbM

˜ ) ------------ --------- --H j j j j j j H b1a1 Hb1a2 · · · Hb1aL Hb1b1 Hb1b2 · · · Hb1bM j j j j j j ··· H ··· H H H H H b2a1

b2a2

b2aL

b2b1

b2b2

b2bM

··· ··· ··· ··· ··· ··· ··· ··· j j j j j j HbMa1 HbMa2 · · · HbMaL HbMb1 H bMb2 · · · HbMbM

The isolated polarizabilities of molecules a and b are aL

π j (a) )

bM

∑ ∑ Hj ij

]

i)a1 j)a1 bM

π j (b) )

bM

∑ ∑ Hj ij

i)b1 j)a1

We now proceed to divide the above 3 × (L + M) by 3 × (L + M) T interaction matrix into four blocks: a 3 × L by 3 × L block, a 3 × L by 3 × M block, a 3 × M by 3 × L block, and a 3 × M by 3 × M block. Denoting the inverse of the 3 × L by 3 × L matrix and the 3 × M by 3 × M matrix as F˜ ˜ and G

F˜ )

[

[

j F a1a1 j F a2a1

··· j F aLa1

j F a1a2 · · · j Fa2a2 · · · ··· ··· j F aLa2 · · ·

j j G G b1b1 b1b2 j j G G ˜ ) b2b1 b2b2 G ··· ··· j j G G bMb1

bMb2

···

j F a1aL j F a2aL

··· j F aLaL

]

j G b1bM j G

··· b2bM ··· ··· j ··· G bMbM

]

we have the gas-phase polarizability for molecules a and b

j (a) ) R

aL

aL

∑ ∑ Fj ij

i)a1 j)a1

j (b) ) R

bM

bM

∑ ∑ Gj ij

i)b1 j)b1

Thus, the 3 × (L + M) by 3 × (L + M) atomic T interaction matrix is replaced by the following 3 × 2 by 3 × 2 molecular interaction matrix

Tãb )

[

j -1(a) R

(jI - Rj-1(a)πj (a)) · πj -1(b)

j -1(b)π j (b)) · π j -1(a) (Ij - R

j -1(b) R

]

Obviously, the product of the above 3 × 2 by 3 × 2 matrix and the isolated 3 × 2 by 3 polarizability matrix yields a 3 × 2 by 3 identity matrix

Tãb ·

[ ] [] π j (a)I

π j (b)I

)

jI

jI

In this way, we find corresponding 3 × 3 matrices for each of the four blocks in the 3 × (L + M) by 3 × (L + M) atomic interaction matrix. This ensures that the assignment of elements in the approximated interaction matrix Tãb is exact for a twomolecular-body system! For a many-molecules system, we carry out the AMPT approximation and replace all blocks (the 3 × 5 by 3 × 5 matrices corresponding to anion-anion interactions, the 3 × 22 by 3 × 22 matrices corresponding to cation-cation interactions, and the 3 × 22 by 3 × 5 cation-anion interaction matrices) representing atomic two-body interactions with corresponding 3 by 3 two-body molecular interaction matrices. The full atomic interaction matrix of 3 × 64 × (22 + 5) ) 5184 by 5184 elements is replaced by a molecular interaction matrix of size 3 × 128 by 3 × 128. From our tests on random system configurations, the inverse method based on the AMPT approximation is about 6 times more efficient than the full inversion method and is about 1.5 times faster than the full selfconsistent method. Unlike the original self-consistent and full inversion methods, the AMPT algorithm does not require storage of a huge matrix and, in principle, can be applied to much larger systems. We have made no attempt at using the AMPT method to solve a self-consistent molecular polarizability problem but expect this approach to be very fast since the number of selfconsistent equations that need to be solved is greatly reduced. Figure 2A and B shows results for the calculation of C2(t) and its negative derivative D(t) under the AMPT diagonalization method compared to the self-consistent method based on atomic


Hu et al.

polarizabilities. The overall agreement is very good, except that the peak of D(t) calculated from the AMPT approximation is slightly lower than the corresponding value calculated from the full atomistic self-consistent approach. Excellent agreement between the two approaches was also found for the separated components (not shown here). The AMPT, which is much less memory demanding, should therefore be a very attractive option for the calculation of OKE spectra for large systems. 4. Simulation Details We performed molecular dynamics simulation for the [MOEPY+][DCA-] system, which has recently been studied by Shirota and Castner18 using software GROMACS version 3.3.45,46 Our simulation box consists of 64 pairs of ions. Periodic boundary conditions are employed using the Particle Mesh Ewald (PME) method to treat long-range electrostatic interactions.17,47,48 In order to avoid sharp cutoffs, we used a switching function at 10.0 Å for the Lennard-Jones interactions. The electrostatic model was obtained by fitting to zero-field density functional theory (DFT) calculations using a large basis set. Stretch, bend, torsion, and van der Waals parameters were taken from the OPLS-AA force field.49 Initial configurations for the gas-phase ions were first minimized by quantum computations at the Hartree-Fock level (6-31G** basis set). Values of the fixed charges were determined by a least-squares fit to the electrostatic potential (ESP), calculated by DFT with the B3LYP method50,51 and the ccpVTZ(-f) basis set. Atomic polarizabilities were then generated using the methodology developed by Stern et al.41 where a dipole moment is placed on each of the atoms including hydrogen. One then applies a series of external electrostatic perturbations in the form of dipolar probes consisting of two opposite charges (0.78e-) with a separation of 0.58 Å. For each perturbation, we computed the change in the ESP at a set of grid points, as given by ab initio calculations using DFT with the B3LYP50,51 method and cc-pVTZ(-f) basis set, on the perturbed and unperturbed conformations. The atomic polarizablities were chosen to minimize the mean-square deviation of the changes in the ESP between the quantum mechanical calculation and our model. All ab initio calculations were performed using the Jaguar program.52 The resulting partial charges and polarizabilities for each atom are listed as Tables 4 and 5 in the Appendix. Using the developments described in section 3, we find spherical polarizabilities (R) for [MOEPY+] and [DCA-] to be 13.36 and 5.4 Å3, respectively. These values are quite large, indicating that the [MOEPY+][DCA-] system may have strong dipole-induced-dipole interactions. In the model that we use for atomic polarizabilities,41 pairs of atoms with either stretching (1-2) or bending (1-3) interactions do not contribute to the total dipole-induced-dipole interactions. This model provides advantages over previous formulations where interactions at short distances were not represented accurately.39 Several other models have been developed to correct for this deficiency (Thole model, Stopdi model, and anisotropic atom point dipole model).53–59 A comparison between these may be of intererest in future studies. Because of the slow heterogeneous17,20,48 nature of ILs, sets of very long simulations are required to fully converge the polarizability-polarizability- time correlation functions at all relevant times. It is because of this that we devised a logarithmic storage scheme in which shorter trajectories with very frequent information recording are used to compute the short part of the polarizability-polarizability correlation function, while intermediate and long time parts are obtained from much longer

Figure 3. Comparison of radial distribution functions in polarizable/ nonpolarizable simulations. Atom N3 is the nitrogen in [MOEPY+], while atom N6 is the center nitrogen atom in [DCA-]. Polarizable simulations were carried out using the SIM package developed in Prof. Berne’s Group.61 Atomic polarizabilities are given in the Appendix. The temperature is T ) 300 K.

trajectories in which information is stored less frequently. This allowed us to capture features ranging from intramolecular vibrations to slow long time scale rotational diffusion. Five independent trajectories were run at T ) 400 K for a duration of 60 ns using a time step of 1 fs. This long simulation time is necessary because room-temperature ionic liquids are dynamically heterogeneous, and otherwise, the calculation of collective polarizability would be extensively biased by initial configurations.20,23 Data for analysis were recorded every 0.1 ps. In order to capture the subpicosecond and femtosecond dynamics, we also collected data every 0.01 and 0.001 ps in each case from 300 simulations. The initial configurations for each of these simulations was obtained from coordinates collected every 1 ns from each of our 5 independent 60 ns trajectories. The duration of these 300 runs was 10 and 2 ps for the case of data collection at 0.01 and 0.001 ps, respectively. In order to study the effect of temperature on the OKE spectrum, we also ran simulations at T ) 500 and 600 K. Averaging is easier as one raises the temperature; therefore, at 500 and 600 K, only 20 and 10 ns trajectories were necessary to reasonably converge the decay of our correlation functions. Correspondingly, 100 and 50 simulations were also run to record data every 0.01 and 0.001 ps at 500 and 600 K, respectively. Since our simulations are very long, we used a Berendsen temperature coupling scheme every 1 ns to prevent minor deviations in temperature.60 Since the rescaling of temperatures occurs so infrequently, we expect the dynamics to be almost identical to what would be obtained in the NVE ensemble. It should be emphasized that the simulations used in order to compute the polarizability-polarizability correlation function do not utilize the polarizable force field. The cost of running 60 ns with a polarizable force field for which computer code is not efficiently made parallel is prohibitive. This is common practice in most if not all previously published articles that compute Raman and OKE spectra. The polarizability is used at a later stage during the analysis of the trajectory and in the computation of the correlation function. Since we have produced a polarizable force field, we carried out some test molecular dynamics runs. Structural properties of the liquid are almost identical when a polarizable force field is used (see Figure 3), but the dynamics is somewhat faster in the case of the polarizable force field. 5. Results and Discussion Before we describe the results derived from our calculation of the OKE transient signal and corresponding OKE spectrum,



TABLE 1: Parameters Obtained by Fitting Rotational Correlation Functions of Different Body-Fixed Cationic Polarizability Eigenvectors

can be observed. These are consistent with the apparent two minima in the velocity autocorrelation function. At high temperatures, one of the minima is only incipient, while at lower temperature, its presence is clear. These two features are also present in our OKE spectral densities (see Figure 5) in which the shoulder becomes a peak as the temperature is lowered. At room temperature, the experimental spectrum shows that the peak at around 7 THz has turned into a shoulder, while the shoulder at around 20 THz becomes the main spectral peak at low frequencies. Judging by the overall shift of the power spectrum, it is clear that intermolecular motion appears at lower frequencies when the temperature is increased. This effect is caused by a decrease in Coulombic interaction strength due to the expansion of the system at higher temperatures. The same sort of phenomenon can be observed in the experimental OKE signal when comparing an ionic liquid with its neutral isoelectronic counterpart.18 At any given temperature, the lowfrequency OKE spectrum appears at higher frequency in the case of the IL due to stronger intermolecular interactions. Different from the power spectra, the OKE signal is weighted by the polarization of different atoms, and it contains more dynamical information about inter- and intramolecular interactions. The experimental and computationally obtained temperature-dependent OKE spectra in the time and frequency domains are shown in Figures 5 and 6, respectively. In Figure 5, we compare the experimentally obtained OKE signal with that derived from our computer simulations at 400, 500, and 600 K. Both the experiment and our calculations clearly show a very sharp peak appearing at less than 0.1 ps-1, corresponding to the long time decay of reorientational motion. The experimental peaks above 30 ps-1 (50 and 110 ps-1 in the case of our simulations) correspond to intramolecular vibrations assigned by Shirota and Castner to be a bending mode of the anion and a ring deformation mode of the cation, respectively. It is clear that our force field shows intramolecular vibrations at frequencies that are not exactly the same as those in the experiment. The anionic bending frequency appears to be about 15 THz higher in simulations, while the ring deformation mode is about 10 THz lower. In our simulations, the intramolecular frequencies appear to be independent of the temperature, but the intensity of these peaks is lower at higher temperature. This is in contrast to the intensity of the lower-frequency part of the OKE signal corresponding to intermolecular motion, which increases in intensity as the temperature is raised. The broad distribution of frequencies from several ps-1 to around 30 ps-1 is due to the complex intermolecular dynamics of the IL. The shift of these low-frequency peaks is consistent with what is observed from our power spectra obtained from

vector

temperature (K)

τ1 (ps)

τ2 (ps)

1

400 500 600 400 500 600 400 500 600

55.6 ( 0.7 15.4 ( 0.2 27.4 ( 0.5 14.11 ( 0.03 2.30 ( 0.02 0.85 ( 0.01 15.62 ( 0.03 2.14 ( 0.02 0.795 ( 0.008

783.1 ( 0.1 98.5 ( 0.1 44.8 ( 0.6 109.9 ( 0.2 10.87 ( 0.04 4.06 ( 0.02 107.4 ( 0.1 10.16 ( 0.03 3.88 ( 0.01

2 3

TABLE 2: Parameters Obtained by Fitting Rotational Correlation Functions of Different Body-Fixed Anionic Polarizability Eigenvectors vector

temperature (K)

τ1 (ps)

τ2 (ps)

1

400 500 600 400 500 600 400 500 600

37.2 ( 0.1 4.00 ( 0.03 1.47 ( 0.02 4.87 ( 0.03 0.808 ( 0.004 0.319 ( 0.004 3.52 ( 0.03 0.768 ( 0.004 0.301 ( 0.003

226.09 ( 0.05 29.09 ( 0.02 11.24 ( 0.01 65.2 ( 0.3 6.57 ( 0.04 1.965 ( 0.002 51.2 ( 0.2 5.51 ( 0.03 1.64 ( 0.02

2 3

TABLE 3: Parameters Obtained by Fitting C2(t) temperature (K)

τ1 (ps)

τ2 (ps)

400 500 600

42 ( 1 0.92 ( 0.03 0.62 ( 0.02

359 ( 2 32.3 ( 0.3 12.52 ( 0.1

it is instructive to analyze the velocity autocorrelation function and its Fourier transform19 (the power spectrum of the liquid) which to some extent will share the characteristics of the OKE spectrum. We computed velocity autocorrelation functions N

VAC(t) )

∑

1 〈ν (0) · νi(t) 〉 N i)1 i

where the index i runs over all cations and anions and Vi(t) is their corresponding velocity of the center of mass. The calculated VACs and corresponding power spectra are displayed in Figure 4. The Fourier transformed spectrum (Figure 4B) at high temperature appears to show only one peak at about 5 THz (1 ps-1 ) 1 THz), while as the temperature is lowered, this peak shifts to higher THz, and a shoulder at around 20 THz

Figure 4. Velocity autocorrelation functions (A) and power spectra (B).


Hu et al.

TABLE 4: Coordinates and Polarizabilities of [MOEPY+]a atom

charge (e)

x (Å)

y (Å)

z (Å)

R (Å3)

H1 H2 H3 H4 H5 N6 C7 C8 C9 C10 C11 C12 C13 C14 H15 H16 H17 H18 H19 H20 H21 O22

0.09287854 0.19239150 0.27792361 0.27792361 0.19239150 0.61503816 -0.09464031 0.15873689 -0.09464031 -0.36933925 -0.36933925 -0.13478131 0.00167019 -0.29082167 0.07905346 0.07905346 0.05363432 0.05363432 0.14694203 0.14694203 0.14694203 -0.16159354

4.24938 3.02142 0.67940 0.73368 3.07810 0.64652 2.58598 3.23154 2.55437 1.25823 1.28909 -0.76629 -1.74655 -4.08012 -0.91010 -0.93099 -1.56939 -1.58655 -5.00644 -4.05773 -4.03547 -3.01373

0.31338 2.22325 1.84972 -2.10707 -1.89681 -0.13539 -1.04349 0.18742 1.25326 1.06596 -1.17942 -0.29325 0.05573 0.13384 -1.32578 0.35744 -0.59775 1.09652 -0.04675 1.17690 -0.53295 -0.13368

-0.83530 0.23835 1.03476 -0.21623 -1.06324 0.42806 -0.61148 -0.48113 0.11422 0.56121 -0.14821 0.86639 -0.26108 -0.61568 1.18747 1.72632 -1.13351 -0.59398 -0.07006 -0.95860 -1.48721 0.29386

0.01000 0.46636 0.65773 0.65773 0.46636 2.23174 0.27040 2.04699 0.27040 0.28044 0.28044 0.45371 1.14150 0.57427 0.30893 0.30893 0.38490 0.38490 0.38490 0.38490 0.38490 0.56943

a Polarizabilities were fit according to ref 41. Because a point dipole description is only accurate for larger distances, 1,2- and 1,3-interactions have been omitted in calculating the gas-phase polarizabilities.

TABLE 5: Coordinates and Polarizabilities of [DCA-] atom

charge (e)

x (Å)

y (Å)

z (Å)

R (Å3)

C1 C2 N3 N4 N5

0.50873113 0.50873113 -0.60193173 -0.70776526 -0.70776526

0.00000 0.00000 0.00000 0.00000 0.00000

1.15285 -1.15285 0.00000 2.24667 -2.24667

0.05355 0.05355 0.68568 -0.38873 -0.38873

1.30551 1.30551 0.75419 1.00994 1.00994

the calculation of VAC functions shown in Figure 4 and is caused by the decrease in intermolecular interaction strengths. This has previously been observed in another ionic liquid (1pentyl-3-methylimidazolium bis(trifluoromethanesulfonyl)imide) by Quitevis and co-workers.11 Because of the extreme computational difficulties described in the Introduction section, we have not attempted to obtain the OKE signal at 300 K; nonetheless, our spectra at 400 K qualitatively agree with the experimental data at room temperature and should be appropriate to provide a reasonable explanation to the dynamics explored by OKE spectroscopy in this ionic system. At short times (less than 100 fs), the time domain spectra in Figure 6A resembles the pulse of the laser. Intramolecular vibrational motion contributes to the subpicosecond oscillations

Figure 5. Frequency domain OKE signals at different temperatures (see eq 2). Experimental data at T ) 300 K were provided to us by Prof. Shirota and Prof. Castner (see Figure 3a of ref 18). Units: to convert from ps-1 to cm-1, the units used in ref 34, we have ν/ps-1 × 100 ) ν/cm-1 × 6π.

at all temperatures. It is clear from Figure 6A that these intramolecular frequencies are slightly different in our simulations when compared to those experimentally observed. On this short time scale, we see that the signal at low temperature decays faster than that at high temperature, an effect that is reversed at long times. This can be clearly seen in Figure 6B. The calculated spectra at three different temperatures (400, 500, and 600 K) along with the experiment data at 300 K obtained by Shirota and Castner18 show that the long time decay of the OKE signal is strongly temperature dependent. At lower temperatures (300 and 400 K), the tails of the time domain signals do not decay to zero even after several hundred picoseconds, while at higher temperatures (500 and 600 K), the signal approaches zero at times less than 100 ps. As Cang and co-workers explain,8 to a first order approximation, the OKE transient signal should be almost identical to a rotational correlation function. This statement has been shown to be correct for many traditional polar solvents such as methanol and acetonitrile37 where, after initial collision-induced phenomena occur on a subpicosecond time scale, rotations dominate the change in polarizability as observed from a laboratory frame of reference. As we will show next, the time scale on which this is true for ILs at room temperature is very different from what is observed in conventional dipolar solvents. So-called “collision-induced” or nonrotational phenomena are quite important for the first few hundred picoseconds at 400 K and probably longer at 300 K, reinforcing the notion that cageing and, in general, nondiffusive dynamics persist in these systems for long times when compared to conventional solvents. This is in contrast to the case of



Figure 6. OKE signal at different temperatures; see eq 1. (A) Initial short time window. (B) Long time window. The experimental data were provided to us by Prof. Shirota and Prof. Castner; see Figure 2 of ref 18.

Figure 7. Time correlation functions C2(t) and their molecular and interaction-induced components; see eq 4. All of the time correlations are normalized by C2(t) at the corresponding temperature; (A) T ) 400 K, (B) T ) 500 K, (C) T ) 600 K.

conventional dipolar solvents far from their glass transition temperatures, where these are only important on a subpicosecond regime. 5.1. Analysis of the Time-Dependent OKE Signal and Comparison to Rotational Time Correlation Functions. Because the fs pulse used in OKE spectroscopy is broad in the frequency domain but sharp in the time domain, the microscopic origin of the OKE signal has different molecular and collective origins. In light of recent experiments, it is interesting to study the different contributions to the anisotropic polarizability relaxation coming from rotational and librational dynamics as well as interaction-induced effects by decomposing the total collective polarizability into several components, as described in detail in section 2. The results of the decomposition of the anisotropic part of the system’s polarizability according to eqs 4 and 5 are shown in Figures 7 and 8. At T ) 400 K, the decomposition of C2(t) into a pure molecular component and an interaction-induced component (Figure 7A) shows that the

interaction-induced component makes substantial contribution to the total anisotropic collective polarizability. This is also true at the higher temperatures studied. The cross term CMI is 2 negative, and its absolute value almost equals the sum of CM 2 (t) and CI2. The correlation function C2(t) is much smaller in I magnitude than either of the terms CM 2 (t) and C2. This is because the interaction-induced polarization nearly cancels the effect of the polarization change caused by molecular motions. This partitioning of the polarizability correlation function gives us interesting information, indicating that interplay between two large and opposite in sign contributions are shaping the OKE response, but unfortunately, as explained in section 2, this partitioning does not provide much physical intuition in terms of what is really happening in the liquid. In contrast, it is more instructive to partition the anisotropic polarizability correlation function into (1) a purely rotational contribution stemming both from the molecular and induced parts of the polarizability and (2) a nonrotational contribution.


Hu et al.

Figure 8. Time correlation functions C2(t) and their reorientational and collisional-induced components; see eq 5. All of the time correlations are normalized by C2(t) at the corresponding temperature; (A) T ) 400 K, (B) T ) 500 K, (C) T ) 600 K.

The results of analyzing the polarizability correlation function in this way are shown in Figure 8. The most striking finding is that at 400 K, the full C2(t) resembles CR2 (t) only after 250 ps! This is in contrast to most conventional solvents in which rotational dynamics dominate the OKE signal after intramolecular time scales and picosecond collisional events occur. At room temperature, rotational phenomena should dominate the anisotropic polarizability fluctuations only on an even much longer time scale! This implies that most OKE experiments on ionic liquids recorded with pump-probe delays that are subnanosecond in duration are only in part capturing rotational dynamics, and a big part of what is observed is nondiffusive cageing dynamics. It is interesting to point out that the effect of cancellation between molecular and induced components of the anisotropic polarizability can be readily seen in the projected decomposition shown in Figure 8A. The reorientational contribution after projection, CR2 (t), is nearly 10 times smaller than that of the molecular contribution before projection, CM 2 (t), in Figure 7A. This strong interaction-induced effect is similar to previous findings in acetonitrile by Ladanyi and Liang.37 As shown in Figure 8B and C, at a higher temperatures, the decay of the total anisotropic collective polarizability-time correlation function, C2(t), follows that of CR2 (t) from much shorter times, consistent with the dynamics being faster and the systems being less dynamically heterogeneous.20,23 Due to the binary nature of ionic liquids, it is important to explore the respective contribution of cations and anions to the overall signal. The corresponding contributions are obtained from eq 6. Figure 9 shows the time correlation functions C2(t) and their cationic, anionic, and cross components at the three temperatures investigated. One interesting observation is that the cation-anion cross-correlation term has opposite sign to the total anisotropic collective polarizability time correlation function C2(t) at all times and at all of the temperatures investigated.

Because the OKE signal is directly proportional to the ˙ 2(t) with the laser pulse, we also plot in Figure convolution of -C 10 the negative first derivative of C2(t), D(t) ) -C˙2(t), at the three temperatures studied. Interestingly, the decay of D(t) resembles that of the cationic contribution D+(t) at all times larger than 200 fs. Both the anionic and cross term contributions are nearly zero after 200 fs. The anionic and cross term correlation functions contribute about 20% at very short times (less than 100 fs). One would therefore expect the cationic contribution to dominate the decay of the OKE signal at all temperatures, as is shown in Figure 11. Because the gas-phase polarizability of [MOEPY+] (13.36 Å3) is much larger than that of the anion [DCA-] (5.4 Å3), it is not surprising that the response of the cation dominates the overall OKE signal. This suggests that a choice of cations and anions with different isolated polarizabilities might control their relative contribution to the dynamical response of the liquid. In order to independently establish the time scales for rotational dynamics in this liquid and compare them at different temperatures with the decay of C2(t), we have computed rotational correlation functions with molecular frame directions along the corresponding eigenvectors of the gas-phase polarizability tensor of cations and anions. All rotational correlation functions for cations and anions as well as C2(t) at temperatures of 400 K and above are well fitted by a sum of two exponentials. Tables 1 and 2 show the parameters obtained from fitting rotational correlation functions defined as 〈vectori(0) · vectori(t)〉, where vectori(t) corresponds to each of the three ionic polarizability eigenvectors shown in Figure 12. It is clear from Table 3 that at all temperatures, the decay of C2(t) at long times appears to be faster than the corresponding decay of the slowest cationic rotational correlation function but faster than the other two rotational components of the cation.



Figure 9. Time correlation functions C2(t) and their components; see eq 6. (A) T ) 400 K; (B) T ) 500 K; (C) T ) 600 K.

Figure 10. Negative first derivative of the time correlation functions C2(t) and its components (see eq 7). D stands for -C˙2; (A) T ) 400 K, (B) T ) 500 K, (C) T ) 600 K.

6. Conclusions We have performed molecular dynamics simulation to study the OKE response of the ionic system [MOEPY+][DCA-]. We have also developed a general algorithm to convert atomic polarizabilities into approximate molecular polarizabilities in order to facilitate the computation of the total collective polarizability based on the DID model. Our computed spectra show all of the features in the OKE spectrum experimentally obtained by Shirota and Castner. Our

analysis of the corresponding time correlation functions of the anisotropic collective polarizability shows that interactioninduced or “nonrotational” effects play an important role in the total collective polarization for times as long as 250 ps at 400 K, and we predict this effect to be much more pronounced at lower temperatures. This indicates that not only rotational dynamics but also cageing dynamics is fundamentally important at all times, consistent with the pump-probe delays used by Shirota and Castner. At 300 K, longer delays should be used to


Hu et al.

Figure 11. Total OKE signal and cationic, anionic, and cross contributions; (A) T ) 400 K, (B) T ) 500 K, (C) T ) 600 K.

CAREER#0547640 awarded to C.J.M. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Appendix Figure 12. Body-fixed vectors of (A) the cation, [MOEPY+], and (B) the anion, [DCA-].

fully observe the decay of the OKE signal due to rotational dynamics. Our OKE spectra show that the signature of intermolecular dynamics shifts to a higher frequency when the temperature is decreased. We find that correlated with the relatively large isolated polarizability of the cation, the OKE signal after initial subpicosecond dynamics is mostly due to [MOEPY+]. The time constant for long time decay of C2(t) is smaller than the slowest cationic rotational component but larger than that of the other two rotational cationic components. Acknowledgment. This fairly large amount of theoretical and computational work requires several acknowledgements to people without whom this work would not have been possible. First, we thank Prof. Shirota and Prof. Castner for kindly sending us their original experimental data. We also would like to thank Dr. Byungchan Kim for useful discussions and great help on building the polarizable and nonpolarizable force fields. This work was started when Dr. Xuhui Huang was still a graduate student at Columbia University. His work was supported by NSF Grant # CHE03-13401 awarded to Prof. Bruce J. Berne at Columbia University. C.J.M. would particularly like to thank Prof. B. J. Berne (under who’s mentoring and great vision he began working on RTILs48) for his generosity, his sound advice, and fruitful discussions. This material is based upon work supported by the National Science Foundation under Grant

The partial charges and polarizablities for each atom are listed in Tables 4 and 5. The method to fit these parameters has been explained in section 4. References and Notes (1) Carmichael, A. J.; Seddon, K. R. J. Phys. Org. Chem. 2000, 13, 591–595. (2) Muldoon, M. J.; Gordon, C. M.; Dunkin, I. R. J. Chem. Soc., Perkin Trans. 2001, 2, 433–435. (3) Aki, S. N. V. K.; Brennecke, J. F.; Samanta, A. Chem. Commun. 2001, 413–414. (4) Karmakar, R.; Samanta, A. J. Phys. Chem. A 2002, 106, 6670– 6675. (5) Karmakar, R.; Samanta, A. J. Phys. Chem. A 2002, 106, 4447– 4452. (6) Hyun, B.-R.; Dzyuba, S.; Bartsch, R.; Quitevis, E. J. Phys. Chem. A 2002, 106, 7579–7585. (7) Antony, J. H.; Mertens, D.; Dolle, A.; Wasserscheid, P.; Carper, W. R. ChemPhysChem 2003, 4, 588–594. (8) Cang, H.; Li, J.; Fayer, M. D. J. Chem. Phys. 2003, 119, 13017– 13023. (9) Giraud, G.; Gordon, C. M.; Dunkin, I. R.; Wynne, K. J. Chem. Phys. 2003, 119, 464–477. (10) Ingram, J. A.; Moog, R. S.; Ito, N.; Biswas, R.; Maroncelli, M. J. Phys. Chem. B 2003, 107, 5926–5932. (11) Rajian, J. R.; Li, S.; Bartsch, R. A.; Quitevis, E. L. Chem. Phys. Lett. 2004, 393, 372–377. (12) Saha, S.; Mandal, P. K.; Samanta, A. Phys. Chem. Chem. Phys. 2004, 6, 3106–3110. (13) Mandal, P.; Sarkar, M.; Samanta, A. J. Phys. Chem. A 2004, 108, 9048–9053. (14) Del Popolo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 1744– 1752. (15) Ito, N.; Arzhantsev, S.; Maroncelli, M. Chem. Phys. Lett. 2004, 396, 83–91.

MD Study of 1-Methoxyethylpyridinium Dicyanoamide (16) Chowdhury, P. K.; Halder, M.; Sanders, L.; Calhoun, T.; Anderson, J. L.; Armstrong, D. W.; Song, X.; Petrich, J. W. J. Phys. Chem. B 2004, 108, 10245–10255. (17) Margulis, C. J. Mol. Phys. 2004, 102, 829–838. (18) Shirota, H.; Castner, E. J. Phys. Chem. A 2005, 109, 9388–9392. (19) Urahata, S. M.; Ribeiro, M. C. C. J. Chem. Phys. 2005, 122, 024511/ 1–024511/9. (20) Hu, Z.; Margulis, C. J. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 831–836. (21) Arzhantsev, S.; Jin, H.; Ito, N.; Maroncelli, M. Chem. Phys. Lett. 2006, 417, 524–529. (22) Hu, Z.; Margulis, C. J. Phys. Chem. B 2006, 110, 11025–11028. (23) Hu, Z.; Margulis, C. J. Acc. Chem. Res. 2007, 40, 1097–1105. (24) Xiao, D.; Rajian, J.; Cady, A.; Li, S.; Bartsch, R.; Quitevis, E. J. Phys. Chem. B 2007, 111, 4669–4677. (25) Maginn, E. J. Acc. Chem. Res. 2007, 40, 1200–1207. (26) Shim, Y.; Jeong, D.; Manjari, S.; Choi, M. Y.; Kim, H. J. Acc. Chem. Res. 2007, 40, 1130–1137. (27) Bhargava, B. L.; Balasubramanian, S. Chem. Phys. Lett. 2006, 417, 486–491. (28) Iwata, K.; Okajima, H.; Saha, S.; Hamaguchi, H. Acc. Chem. Res. 2007, 40, 1174–1181. (29) Urahata, S.; Ribeiro, M. J. Chem. Phys. 2006, 124, 074513/1– 074513/8. (30) Castner, E. W.; Wishart, J. F.; Shirota, H. Acc. Chem. Res. 2007, 40, 1217–1227. (31) Hunt, N. T.; Jaye, A. A.; Meech, S. R. Phys. Chem. Chem. Phys. 2007, 9, 2167–2180. (32) Ozawa, R.; Hayashi, S.; Saha, S.; Kobayashi, A.; Hamaguchi, H. Chem. Lett. 2003, 32, 948–949. (33) Saha, S.; Hayashi, S.; Kobayashi, A.; Hamaguchi, H. Chem. Lett. 2003, 32, 740–741. (34) Shirota, H.; Castner, E. J. Phys. Chem. B 2005, 109, 21576–21585. (35) Frenkel, D.; McTague, J. P. J. Chem. Phys. 1980, 72, 2801–2818. (36) Geiger, L. C.; Ladanyi, B. M. Chem. Phys. Lett. 1989, 159, 413– 420. (37) Ladanyi, B. M.; Liang, Y. Q. J. Chem. Phys. 1995, 103, 6325– 6372. (38) Ryu, S.; Stratt, R. M. J. Phys. Chem. B 2004, 108, 6782–6795. (39) Applequist, J.; Carl, J. R.; Fung, K. J. Am. Chem. Soc. 1972, 94, 2952–2960. (40) Murry, R. L.; Fourkas, J. T.; Keyes, T. J. Chem. Phys. 1998, 109, 2814–2825.

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7849 (41) Stern, H.; Kaminski, G.; Banks, J.; Zhou, R.; Berne, B.; Friesner, R. J. Phys. Chem. B 1999, 103, 4730–4737. (42) McMorrow, D.; Lotshaw, W. T.; Kenney-Wallace, G. A. IEEE J. Quantum Electron. 1989, 24, 443–454. (43) McMorrow, D.; Lotshaw, W. T. J. Phys. Chem. 1991, 95, 10395– 10406. (44) Ladanyi, B. M. J. Chem. Phys. 1983, 78, 2189–2203. (45) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43–56. (46) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (47) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577–8593. (48) Margulis, C. J.; Stern, H. A.; Berne, B. J. J. Phys. Chem. B 2002, 106, 12017–12021. (49) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (50) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785–789. (51) Becke, A. D. Phys. ReV. A 1988, 38, 3098–3100. (52) Jaguar, version 3.5; Schrodinger, Inc.: Portland, OR, 1998 (http:// www.schrodinger.com). (53) Birge, R. R. J. Chem. Phys. 1980, 72, 5312–5319. (54) Birge, R. R.; Schick, G. A.; Bocian, D. F. J. Chem. Phys. 1983, 79, 2256–2264. (55) Dang, L. X.; Chang, T. M. J. Chem. Phys. 1997, 106, 8149–8159. (56) Elola, M. D.; Ladanyi, B. M. J. Chem. Phys. 2005, 122, 224506/ 1–224506/15. (57) Paolantoni, M.; Ladanyi, B. M. J. Chem. Phys. 2002, 117, 3856– 3873. (58) Thole, B. T. Chem. Phys. 1981, 59, 341–350. (59) van Duijnen, P. T.; Swart, M. J. Phys. Chem. A 1998, 102, 2399– 2407. (60) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (61) Stern, H.; Xu, H.; Harder, E.; Rittner, F.; Pavesse, M.; Berne, B. J. SIM: A Molecular Dynamics Simulation Program. (62) Yan, T.; Li, S.; Jiang, W.; Gao, X.; Xiang, B.; Voth, G. A. J. Phys. Chem. B 2006, 110, 1800–1806. (63) Yan, T.; Burnham, C. J.; Del Pópolo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 11877–11881.

JP800729G

Molecular Dynamics Study of the Temperature-Dependent Optical

Recommend Documents