Guest Dynamics in Endohedrally Doped Fullerenes - American

of up to ∼13 cm-1. Analysis of the metal motions inside the cage shows that they can be chaotic. The origin of such chaotic behavior is discussed. I...
0 downloads 0 Views 50KB Size
15048

J. Phys. Chem. B 2005, 109, 15048-15051

Guest Dynamics in Endohedrally Doped Fullerenes Ewout Kemner and Francesco Zerbetto* Dipartimento di Chimica G. Ciamician, UniVersita` di Bologna, V. F. Selmi 2, 40126 Bologna, Italy ReceiVed: April 29, 2005; In Final Form: June 16, 2005

We use molecular dynamics to calculate the vibrational density of states of several endohedral monometallofullerenes and compare the results with the experimental values. The vibrational patterns depend weakly on the charge transfer from the metal to the fullerene cage, with the largest variations of ∼2 cm-1 observed upon switching off completely the Coulomb interactions, and more strongly on the cage isomer, with variations of up to ∼13 cm-1. Analysis of the metal motions inside the cage shows that they can be chaotic. The origin of such chaotic behavior is discussed.

Introduction Endohedrally doped fullerenes are attracting interest both for their chemically nontrivial architecture and for their potential use in molecule-based devices.1 The incorporation of metal atoms inside a fullerenic cage results in strong atomic confinement and in charge transfer from the metal to the carbon structure, with the formation of a stable and tightly held ion pair. A further consequence of this type of doping is the creation of low-lying electronic states delocalized on the cage, which are important for the production of large optical nonlinearities and photoinduced reactions.2 The electronic properties of endohedral fullerenes are a function of the motions of the encapsulated species. Indeed, such properties are formally obtained by averaging over snapshots of the guest dynamics. During the metal orbit, the charge transfer likely goes the carbon atoms that happen to be the closest. Understanding the guest’s motions is therefore desirable, and several mono-metallofullerenes have been studied before by a variety of techniques.3-12 Infrared, Raman, and inelastic neutron scattering experiments together with simulations were used to probe the vibrational motions of a variety of endohedrally doped fullerenes that ranged from La@C82 to Ca@C82, to Y@C82, to Ce@C82, to Gd@C82, to Tm@C82, and to Eu@C74. Two features of this work are important: on one hand, consensus has emerged that the order of stability of the empty and filled cages may be substantially different and that different dopants appear in different cages (i.e., three stable Tm@C82 isomers11a and four stable Ca@C82 isomers11b are found); on the other hand, it is less evident how this structural variability affects the motions of the dopant inside the cage. Standard computational techniques evaluate atomic motions in the harmonic approximation. This is done through the calculation of the vibrational frequencies, and the associated normal modes, from the Hessian matrix (the matrix of the second derivatives of the energy with respect to the nuclei displacements). The implicit assumption of this approach is that there is a single dominant equilibrium position about which all the motions take place. In the case of the endohedrally doped fullerenes, this is hardly satisfactory since the entire inner space of the cage is usually explored by the dopant and a large number of easily interconverted equivalent minima exists. * Corresponding author. E-mail: [email protected].

The exploration of the motions of the endohedral dopant should therefore start with molecular dynamics so that all the minima and interconverting structures are properly explored. Here, we investigate, with a molecular dynamics approach, several endohedral mono-metallofullerenes. From the dynamics, we calculate their vibrational patterns and compare the results with those found in the literature. This enables us to assess the accuracy of the model. We further examine closely the motions of the endohedral atom and find that its orbit can vary greatly in nature. Computational Background To investigate the dynamics of endohedral fullerenes, we used molecular dynamics and Fourier transformed the atom velocity autocorrelation functions to obtain the vibrational patterns13 of M@C82 (M ) La, Y, Gd, Ce, and Tm) and Eu@C74. The vibrations with a low frequency have a strong contribution from the metal motion and are compared with those found in the literature. All calculations were performed on individual endohedrally doped cages with the Tinker program14 and the MM3 force field.15 This computational approach has found many applications in our laboratory.16 We start the dynamics from optimized structures, using, initially, the C82(9) isomer. The time step was 1 fs, and velocities were collected every 5 fs during 500 ps in the NVE ensemble after equilibration at room temperature. Unless specified, the electronic structure of endohedral fullerenes is simply described as M2+/3+@C2n2-/3-, where the charge on the cage is spread evenly on all the carbon atoms. In practice, for atom k, the density of vibrational states at frequency ω is given by

gk(ω) ) C

∫ exp(iωt)〈νk(0)νk(t)〉dt

(1)

where C is a normalization constant such that ∫gk(ω)dω ) 1, and 〈νk(0)νk(t)〉 is the velocity autocorrelation function of the atomic motions; the brackets indicate ensemble averaging. This quantity oscillates in time as the thermal energy flows between the degrees of freedom of the system (i.e., the vibrational levels). At each frequency, ω, one can extract the contribution to 〈νk(0)νk(t)〉 of the vibrational states by a projection over the corresponding oscillating wave of the form sin ωt (Fourier transform). A similar approach has been used before to calculate

10.1021/jp0522213 CCC: $30.25 © 2005 American Chemical Society Published on Web 07/14/2005

Guest Dynamics in Endohedrally Doped Fullerenes

J. Phys. Chem. B, Vol. 109, No. 31, 2005 15049

the inelastic neutron scattering of interlocked molecules17 and exohedrally doped fullerenes.18 To quantify the chaotic behavior of the orbit of the endohedral atom, we calculated the Lyapunov exponents. They are a simple way of discriminating between types of orbits and measuring chaos. One can consider two initial points, xi and xi + ∆x0, for the position of the endohedral atom. The atomic orbits from either point can be very similar or diverge substantially, even when ∆x0 is small. The Lyapunov exponent is defined as

λ)

1 |∆x(xi,t)| ln tf∞|∆x0|f0 t |∆x0| lim

(2)

If λ < 0, the orbit is stable. Negative Lyapunov exponents are characteristic of dissipative or nonconservative systems. The more negative the exponent, the greater the stability. If λ ) 0, the orbit is a neutral fixed point, and the system is in a steady state and conservative. If λ > 0, the orbit is unstable and chaotic. Nearby points, no matter how close, diverge. All neighborhoods in the phase space are eventually visited. This does not preclude the emergence of a pattern in the orbit. A physical example is the Brownian motion. Although the system is deterministic, there is no order to the orbit.

TABLE 1: Experimental and Calculated Vibrations of Endohedral Fullerenesa

system

experimental wavenumbers (cm-1)

Y@C82(9)

54b 85b 150b 183b

70b 128b 150b 180b

La@C82(9)

45-50b

55b, 50c, 27d, 30d 110b, 80c 150b, 159d 180b, 160-180c

163b Ce@C82(9) 162b Gd@C82(9) 151e, 155b Tm@C82(9)

42f

Results and Discussion The investigation of the motions of endohedrally doped fullerenes has been carried out experimentally and theoretically.6 Lebedkin3 et al. studied M@C82 using a variety of techniques: with Raman spectroscopy M was La, Y, Ce, and Gd; with farinfrared spectroscopy, M was La, Y, and Ce; with inelastic neutron scattering (INS), M was La and Y. Krause et al. used Raman and infrared spectroscopy to obtain the spectra of Tm@C824,5 and [email protected] Kuran et al. studied Eu@C74 with infrared and Raman spectroscopy.6 Fullerenes can have different isomers that satisfy the isolated pentagon rule. For the present systems, C74 has only one isomer, while there are nine isomers of C82. Following Fowler and Manolopoulos,7 they are C82(1), C82(3), and C82(5) with C2 symmetry; C82(2), C82(4), and C82(6) with Cs symmetry; C82(7) and C82(8) with C3V symmetry; and finally, C82(9) with C2V symmetry. 13C NMR experiments show that the most commonly observed structure for C82 has C2 symmetry.8 Calculations show that C82(3) is the most stable isomer.9 The endohedral doping of fullerenes may change the relative order of stability of the isomers.10 Calculations show that C82(9) is the most stable isomer for C822- and C823-,9 while it is known experimentally that M@C82 gives several isomers.11 Considering the relative abundance of data, we start with the assessment of the accuracy of the model. Table 1 compares the results of calculations and experiments. A typical example of the calculated vibrational density of states is shown in Figure 1 for La@C82(9). The intensities of the bands reflect the population of the motion at a given wavenumber during the molecular dynamics simulation. Analysis of the velocity autocorrelation function shows that the three lower frequency bands correspond to cage vibrations, while the three higher frequency bands arise from the coupling of the motions of the metal atom and the cage atoms. As expected, the pattern of frequencies shifts towards lower values when the metal mass increases. There appears to be a systematic overestimate of the lowest frequency (actually wavenumber) often located experimentally at ∼50 cm-1, while the higher frequency metal motions are predicted rather accurately. More in detail for the systems investigated here, one can notice that (a) Y@C82: far-infrared absorption experiments find

calculated wavenumbers in previous work (cm-1)

116d, 116-120f Eu@C74 a

c

123g 225g

calculated wavenumbers in the present work (cm-1) 82, 85 96 169 180, 182 189 76,78 89 155 167, 178 191 74,76 87 151 163, 175 190 70,72 82 142 153, 164 190 67,68 78 133 146, 156 192 101, 103, 114 203

In bold, the calculated value of the first cage vibration. Ref 19. d Ref 20. e Ref 5. f Ref 4. g Ref 6.

b

Ref 3.

Figure 1. Calculated vibrational density of states for La@C82(9) at 300 K.

bands at 54 and 180.5 cm-1. Raman identifies a band at 183 cm-1. INS measurements locate bands at 85, 150, and 180 cm-1. The present calculations give bands at 82 and 85, at 96, at 169, and at 180 and 182 cm-1, with cage motions starting at 189 cm-1. The bands at 82, 85, and 96 cm-1 involve motions of the metal atom only. The bands at 169, 180, and 182 cm-1 arise from coupling of the motions of the metal atom and the fullerene cage. Because of the overestimate mentioned previously, the bands calculated at 82 and 85 cm-1 should correspond to the experimental band at 54 cm-1. The metal vibration calculated at 96 cm-1 corresponds to the band measured at 85 cm-1. The band at 169 cm-1 corresponds to the experimental band at 150 cm-1, and the two calculated bands at 180 and 182 cm-1 are assigned to the measured band around 183 cm-1. Previously, with a classical lattice-dynamical model, but for the C82(3) isomer and not the C82(9) isomer, Lebedkin et al. calculated vibrations around 70, 128, 150, and 180 cm-1, where the first two vibrations contain the main metal contributions.3 This agrees with our results: the band of Lebedkin et al. at 70 cm-1 corresponds to our bands at 82 and 85 cm-1, and the band at 128 cm-1 corresponds to our metal band at 96 cm-1. The

15050 J. Phys. Chem. B, Vol. 109, No. 31, 2005 bands at 150 and 180 cm-1 correspond to our bands at 169 and at 180 and 182 cm-1. (b) La@C823: far-infrared absorption experiments show bands around 45-50 and at 161.5 cm-1. Raman scattering and INS find a band at 163 cm-1. Interestingly, the bands found with INS for Y@C82 at 85 and 150 cm-1 do not seem to have a clear counterpart for La@C82. Our calculations give metal vibrations at 76, 78, and 89 cm-1 and metal-cage coupled modes at 155, 167, and 178 cm-1. The internal cage modes start at 191 cm-1. The first two metal vibrations are assigned to the band measured around 45-50 cm-1. The calculated bands at 89 and 155 cm-1 have no experimental counterpart, probably because they are forbidden experimentally by selection rules. Finally, the two calculated bands at 167 and 178 cm-1 correspond to the experimental band around 163 cm-1. Computationally, La@C82 has been studied by several methods. With the classical lattice-dynamical model, also used for Y@C82, Lebedkin et al. determine bands at 55, 110, 150, and 180 cm-1 (for the C82(3) isomer).3 The two lowest frequency bands involve mainly metal atom motions. With ab initio molecular dynamics, Andreoni et al.19 calculated, for the C82(3) isomer, metal vibrations at 50 and 80 cm-1, while at 160180 cm-1 they found motions of La coupled to those of carbon atoms. Note that these results agree reasonably well with the present results: we calculate metal modes at 76, 78, and 89 cm-1 and metal-cage modes at 155, 167, and 178 cm-1. Kobayashi et al.20 calculated the vibrational spectrum of La@C82 using density functional theory and the C82(9) isomer. They found bands at 27, 30, and 150 cm-1. (c) Ce@C82: for this species, the measurements are less exhaustive. Both Raman scattering and far-infrared absorption give a band at 160-162 cm-1.3 This vibration is satisfactorily reproduced by the calculated motions at 163 or 175 cm-1. (d) Gd@C82: Raman spectroscopy locates vibrations at 155 cm-1 (ref 3) and at 151 cm-1.5 We assign this band to the calculated bands at 153 and 164 cm-1. (e) Tm@C82: three stable isomers have been isolated. Depending on the isomer, Raman spectroscopy shows bands at 25 and 42 and at 116 cm-1,4,5 while infrared spectra show bands at 116-120 cm-1. The lowest band, at 25 cm-1, is attributed to a lattice mode. Our calculations give the familiar pattern. The metal vibration modes are present at 67, 68, and 78 cm-1. The bands at 133, 146, and 156 cm-1 arise from metal-cage motion coupling, and the first cage mode is calculated at 192 cm-1. The bands at 67 and 68 cm-1 correspond to the band measured at 42 cm-1. The next two calculated bands, at 78 and 155 cm-1, do not have a Raman or infrared counterpart. The highest experimental band, at 116 cm-1, is attributed to metalcage coupling and is calculated at 146 and 156 cm-1. For this system, the difference between the experimental metal-cage coupling band and our calculations is larger than found previously for the other M@C82. The exceptionally low frequency of this metal excitation has been noted before from an analysis of the experimental band position plotted against the square root of the reciprocal metal mass for various endohedral systems,5 where the band position for Tm@C82 was lower than expected. This was attributed to the different charge of Tm (+2 instead of +3), resulting in a weaker interaction. Our calculations do not suggest that this is the case since they are performed with a charge of +2 for Tm. The reason of the deviation of the experimental value from those of other endohedrally doped fullerenes remains an open question.

Kemner and Zerbetto TABLE 2: Calculated Vibrations, cm-1, for Four La@C82 Isomers (see text) La@C82(1)

La@C82(3)

La@C82(5)

La@C82(9)

63 96 100 159 178

68 86 100 174 181

72 85 85 163, 170, 176, 180 186

76 78 89 155, 167, 178 191

(f) Eu@C74: the Raman spectrum of Eu@C74 shows a bands at 123 cm-1 that are attributed to a metal-C74 cage mode.6 The cage vibrations start at 225 cm-1. Our calculations give frequencies at 101, 103, and 114 cm-1. The first cage mode is found at 203 cm-1. The overestimate of the lowest frequencies makes the model stiffer than reality, a feature that is important to bear in mind in the discussion of the chaotic behavior of the guest motions. Two issues are worth examining because of their relevance to the metal dynamics. They are (i) the influence of the charge transfer and (ii) the effect of different isomeric cages. To investigate the influence of charge transfer on the calculated frequencies, we recalculated the vibrational density of states of La@C82(9) without Coulomb interactions. The new positions (77, 87, 153, 165, 176, and 191 cm-1) are close to those found when the Coulomb interactions are present (7678, 89, 155, 167, and 178 cm-1). Furthermore, we examined the role of a variation of the charge transfer with the dynamics. To calculate the charges of the individual atoms on-the-fly as the metal rattles inside the cage, we used the charge equilibration method of Goddard and Rappe21 that was used in our group for other applications.22-24 No significant change in the positions of the frequencies was observed. In general, the electrostatic interactions must be present to reproduce the off-center position of metal atoms inside the fullerene cages,20 with the metal ions attracted by regions of high electron density in the carbon cage. However, due to the very small charges present on each carbon atom, the Coulomb component of the potential energy surface is rather flat and, as such, has negligible influence on the vibrational spectrum. Indeed, it was shown by Kobayashi et al. that the bottom of the electrostatic potential map of C823- is flat.9 To investigate the influence of the cage isomer on the frequencies, we considered four La@C82 (they are C82(1), C82(3), and C82(5), all with C2 symmetry, and C82(9), with C2V symmetry). The results are shown in Table 2. The lowest frequency cage mode shifts by ∼8 cm-1 for the different isomers. The metal-cage modes vary by up to 13 cm-1 for the first mode, 18 cm-1 for the second mode, and 11 cm-1 for the third. It is concluded that the guest dynamics is more sensitive to variations of the shape of the cage than to the amount of charge transfer. Interestingly, a study of Ne@C60 showed that the dynamics of Ne inside the C60 cage is nonlinear.25 We therefore decided to look at the orbits of the endohedral metals with a different approach. Qualitatively, the motion of the metal can reside on a periodic orbit or can depart from it; the departure is ascribable to the cage vibrations that disturb randomly the motion of the endohedral atom. It was thought that the disturbance could generate chaotic orbits, which are quantifiable through the Lyapunov exponent, λ. This exponent gives the rate of exponential divergence from perturbed initial conditions and is a measure of chaos. For a periodic motion, λ is zero; for a chaotic motion, λ is positive. A negative λ means the system gravitates toward a stable point in phase space.26 Here, we use the method originally proposed by Benettin et al.,27 in which

Guest Dynamics in Endohedrally Doped Fullerenes

J. Phys. Chem. B, Vol. 109, No. 31, 2005 15051

TABLE 3: Largest Lyapunov Exponents λ (ps-1) for the Motion of the Metal Atoms in C82(9) at 300 Ka system

λb

λc

metal

Sc@C82 11.7 0.00046 Sc Y@C82 9.2 0.00039 Y La@C82 8.3 0.00042 La Ce@C82 7.8 0.00036 Ce Gd@C82 7.4 0.00040 Eu Tm@C82 9.8 0.00035 Gd Eu@C74 5.2 0.00036 Tm

mass

radius

44.956 88.905 138.906 139.905 152.921 157.924 168.934

2.61 2.71 2.78 2.74 2.94 2.71 2.67

energy-well depth charge 0.174 0.250 0.312 0.340 0.271 0.365 0.393

+2 +3 +3 +3 +2 +3 +2

a Ion masses (Da), radii (Å), energy-well depths (kcal mol-1), and charges used in the simulations are also given. b Full dynamics. c Dynamics with a rigid cage.

regimes, which are indicated by the two linear fits on the data. Up to 50 K, λ increases slowly with the temperature. From ∼50 K, λ increases sharply, indicating the onset of the real chaotic process in the dynamics. Conclusion We have explored the dynamics of endohedrally doped fullerenes with a model based on molecular dynamics at room temperature. The vibrations due to the metal motions are split into two sets with the lower frequency movements of the guest predicted to be too fast. The cage structure affects the value of the frequency more than the presence of the local charges transferred from the metal to the carbon frame. The charges, however, strongly contribute to the nonlinear, chaotic dynamics of the ions that occurs from above 50 K. Acknowledgment. Support from the European Union Projects HPRN-CT-2002-00177, WONDERFULL, is gratefully acknowledged. F.Z. also thanks MURST project FIRB Carbonio micro e nanostrutturato. References and Notes

Figure 2. Logarithm of the largest Lyapunov exponent, λ, for the La motions in La@C82(9) vs the inverse of temperature; T is in K. The lines are linear fits that indicate the two dynamical regimes that cross over ∼50 K.

the initial configuration is slightly perturbed and its propagation is followed in time. The time step was 0.1 fs, the norm of the initial perturbation was 0.0001 Å, and λ was determined at each time step during 5 ps, after which the values of λ were averaged. The results for different endohedral fullerenes are shown in Table 3. For comparison, we also include the result for Sc@C82. In all cases, the dynamics is nonlinear, although to a varying degree. Qualitatively, one can expect that light atoms with smaller radii are more prone to chaotic vagaries than big, heavy atoms. More strongly bounded (with higher frequency motions) atoms should have smaller λ than more loosely bounded ones. Atoms in a smaller cage are more constrained and should have a smaller λ. Inspection of Table 3 proves that this qualitative reasoning is correct. The lightest atom, Sc, which also has the smallest radius, shallowest energy-well depth, and smallest charge, has the largest λ. Y, La, Ce, and Gd (in order of weight) have a +3 charge and energy-well depths that increase with increasing weight. Eu@C74 has a lower λ because of the confinement in a smaller cage. Tm is not an exception: it is the heaviest atom with the deepest energy well, but its fairly large λ value is due to its lower charge. A test run where Tm was assigned a charge of +3 indeed reduced λ to 7.4, in accordance with the general trend described previously. While Coulomb interactions hardly influence the vibrational pattern, they significantly affect the nature of the orbit. To understand the origin of the chaotic motion, λ values were recalculated with the carbon atoms of the cage fixed in their equilibrium positions; in practice, the carbon atoms only provided the force experienced by the endohedral atom. The results are also given in Table 3, where one can see that all λ values now approach zero. This indicates the existence of a real periodic motion. The slightly positive values are due to the numerical inaccuracy of the algorithms. Naively, one can think that the cage vibrations kick the metal atom and give rise to the nonlinear behavior. Figure 2 shows λ versus 1/T for La@C82(9) at different temperatures. At low temperatures, λ is smallest. In this Arrhenius-type of plot, one can distinguish two dynamics

(1) Bethune, D. S.; Johnson, R. D.; Salem, J. R.; de Vries, M. S.; Yannoni, C. S. Nature 1993, 366, 123. (2) Liu, S.; Sun, S. J. Organomet. Chem. 2000, 599, 74. (3) Lebedkin, S.; Renker, B.; Heid, R.; Schober, H.; Rietschel, H. Appl. Phys. A 1998, 66, 273. (4) Krause, M.; Hulman, M.; Kuzmany, H.; Kuran, P.; Dunsch, L.; Dennis, T. J. S.; Inakuma, M.; Shinohara, H. J. Mol. Struct. 2000, 521, 325. (5) Krause, M.; Kuran, P.; Kirbach, U.; Dunsch, L. Carbon 1999, 37, 113. (6) Kuran, P.; Krause, M.; Bartl, A.; Dunsch, L. Chem. Phys. Lett. 1998, 292, 580. (7) Fowler, P. W.; Manolopoulos, D. E. An atlas of fullerenes; Oxford University Press: Oxford, 1995. (8) Achiba, Y.; Kikuchi, K.; Aihara, Y.; Wakabayashi, T.; Miyake, Y.; Kainosho, M. Mater. Res. Soc. Symp. Proc. 1995, 359, 3. (9) Kobayashi, K.; Nagase, S. Chem. Phys. Lett. 1998, 282, 325. Orlandi, G.; Zerbetto, F.; Fowler, P. W. J. Phys. Chem. 1993, 97, 1357513579. (10) Fowler, P. W.; Zerbetto, F. Chem. Phys. Lett. 1995, 243, 36-41. (11) Kirbach, U.; Dunsch, L. Angew. Chem. 1996, 108, 2518. Xu, Z.; Nakane, T.; Shinohara, H. J. Am. Chem. Soc. 1996, 118, 11309. (12) Andreoni, W.; Curioni, A. Appl. Phys. A 1998, 66, 299. (13) Brawer, S. J. Chem. Phys. 1983, 79, 4539. Goncalves, S.; Bonadeo, H. Phys. ReV. B 1992, 46, 10738. (14) Ponder, J.; Richards, F. J. Comput. Chem. 1987, 8, 1016. Kundrot, C.; Ponder, J.; Richards, F. J. Comput. Chem. 1991, 12, 402. Dudek, M. J.; Ponder, J. J. Comput. Chem. 1995, 16, 791. (15) Allinger, N. L.; Yuh, Y. H.; Lii, J.-H. J. Am. Chem. Soc. 1989, 23, 8551. Lii, J.-H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 23, 8566. Lii, J.-H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 23, 8576. (16) Acocella, A.; Venturini, A.; Zerbetto, F. J. Am. Chem. Soc. 2004, 126, 2362. Gournis, D.; Georgakilas, V.; Karakassides, M. A.; Bakas, T.; Kordatos, K.; Prato, M.; Fanti, M.; Zerbetto, F. J. Am. Chem. Soc. 2004, 126, 8561. Leigh, D. A.; Venturini, A.; Wilson, A. J.; Wong, J. K. Y.; Zerbetto, F. Chem.sEur. J. 2004, 10, 4960. Melle-Franco, M.; Kuzmany, H.; Zerbetto, F. J. Phys. Chem. B 2003, 107, 6986. (17) Bottari, G.; Caciuffo, R.; Fanti, M.; Leigh, D. A.; Parker, S. F.; Zerbetto, F. ChemPhysChem 2002, 3, 1038. Leigh, D. A.; Parker, S. F.; Timpel, D.; Zerbetto, F. J. Chem. Phys. 2001, 114, 5006. (18) Kemner, E.; Zerbetto, F. Chem. Phys. Lett. 2005, 405, 270. (19) Andreoni, W.; Curioni, A. Phys. ReV. Lett. 1996, 77, 834. (20) Kobayashi, K.; Nagase, S. Mol. Phys. 2003, 101, 249. (21) Rappe, A. K.; Goddard , W. A., III. J. Phys. Chem. 1991, 95, 3358. (22) Baxter, R. J.; Teobaldi, G.; Zerbetto, F. Langmuir 2003, 19, 7335. (23) Montalti, M.; Prodi, L.; Zaccheroni, N.; Baxter, R. J.; Teobaldi, G.; Zerbetto, F. Langmuir 2003, 19, 5172. (24) Baxter, R. J.; Rudolf, P.; Teobaldi, G.; Zerbetto, F. ChemPhysChem 2004, 5, 245. (25) Bug, A. L. R.; Wilson, A.; Voth, G. A. J. Phys. Chem. 1992, 96, 7864. (26) Eckmann, J.-P.; Ruelle, D. ReV. Mod. Phys. 1985, 57, 617. (27) Benettin, G.; Galgani, L.; Strelcyn, J.-M. Phys. ReV. A 1976, 14, 2338.