Semiempirical Molecular Dynamics Study of Empty C2(3)-C82

Oct 19, 2009 - Semiempirical Molecular Dynamics Study of Empty C2(3)-C82 Fullerene in Neutral and Charged Forms: Geometrical and Spectroscopic ...
0 downloads 0 Views 3MB Size
19658

J. Phys. Chem. C 2009, 113, 19658–19663

Semiempirical Molecular Dynamics Study of Empty C2(3)-C82 Fullerene in Neutral and Charged Forms: Geometrical and Spectroscopic Characterization Roland Sˇolc,† Vladimı´r Lukesˇ,*,† Michal Ilcˇin,† Peter Rapta,†,‡ Michal Zalibera,†,‡ and Lothar Dunsch*,‡ Institute of Physical Chemistry and Chemical Physics, SloVak UniVersity of Technology, Radlinske´ho 9, SK-81 237 BratislaVa, SloVakia, and Department of Electrochemistry and Conducting Polymers, Leibniz-Institute of Solid State and Materials Research, IFW Dresden, Helmholtzstrasse 20, D-01069 Dresden, Germany ReceiVed: June 29, 2009; ReVised Manuscript ReceiVed: September 26, 2009

A systematic and comparative theoretical calculation based on the Density Functional Theory and semiempirical Austin Model 1 (AM1) and Zerner’s Intermediate Neglect Overlap (ZINDO) methods have been performed on the neutral and charged states of empty C2(3)-C82 fullerene. The relations between the dynamical structure and the total molecular charge were studied using AM1 molecular dynamics (MD). The obtained results show that the shape of the fullerene ball is more-or-less periodically changed around the most probable geometry. These changes were quantified by using the averaged geometrical bond length alternation parameter. On the basis of the evaluated MD geometries and relevant ZINDO transitions, the absorption spectra were simulated. These simulations allowed us to assign individual bands in spectroelectrochemical records where the mixtures of oxidation and reduction products occurred. 1. Introduction Fullerenes represent a unique class of carbon compounds with a regular three-dimensional architecture1,2 the stability of which is controlled by the isolated pentagon rule (IPR).3 While the most abundant and most studied empty fullerenes C60 and C70 exist only in a single isomeric form, nine isomers obeying the IPR are possible for the C82 cage.3 However, until today, only three C82 isomers have been extracted from the fullerene soot in a significant abundance.4 The most stable isomer C82(3) with C2 symmetry was recently isolated in our group and studied in detail by cyclic volammetry and spectroelectrochemistry.5,6 The experimental results of this study, including quantum chemical time-dependent density functional theory (TD-DFT)7 calculations, confirmed the assignment of the structure of the abovementioned isomer. Special efforts were also devoted to the generation and characterization of the charged states from monocation (+1) to tetraanion (-4). Remarkable differences in the ESR 13C satellites pattern of the radicals were found. The largest spread of hyperfine constants (hfc) was observed for the monocation radical. Although the mono- and trianions of C2(3)-C82 have a similar spread of hfc, the distribution pattern is different. To understand the relations between the dynamical structure, total molecular charge, and spectroscopic properties in more detail, the above-mentioned static quantum chemical calculations need to be extended. In this context, the classical molecular dynamics simulations (MD) can significantly contribute to this issue.8 However, each level of MD calculation has its own advantages and disadvantages. The ab initio or DFT levels have a high chemical precision and a few free optional parameters, but the time and size scales of such simulations are extremely * To whom correspondence should be addressed. E-mail: vladimir.lukes@ stuba.sk and [email protected]. † Slovak University of Technology. ‡ Leibniz-Institute of Solid State and Materials Research.

restricted by their computational complexity.9 On the other hand, the force fields methods make a calculation of systems containing a large number of atoms (more than 1000) possible. Unfortunately, the specific effects on electronic state changes are not accounted for. One way of overcoming these computational and theoretical limits is the use of MD based on the semiempirical Hamiltonian. This enables the evaluation of long nanosecond time-scale dynamics in real time, which is a crucial step for the reliable description of vibronic motion or other thermodynamic quantities. The semiempirical methods based on the neglect of diatomic differential overlap (NDDO)10 approximation have been successfully used in quantum chemical studies of various carbon clusters.11-14 Taking the above arguments into account, we performed theoretical studies in connection with the spectroelectrochemical experiments of the C2(3)-C82 fullerene (see Figure 1) in its neutral as well as its charged states. The molecular motions and structural changes during the time evaluation were simulated by the application of the MD, which utilized the semiempirical Austin Model 1 (AM1) Hamiltonian.15 On the basis of the obtained set of MD geometries, the absorption bands were simulated for the neutral as well as charged states. Comparing to the static calculations of excitation energies, these simulations take into account the natural broadening of the absorption bands due to the vibronic motion and therefore the general shape of the simulated spectrum might be closer to the experimental results. Finally, the theoretical results were compared with the reference DFT calculations as well as the experimental results. 2. Computational Details The optimal geometries of the C2(3)-C82 fullerene for the neutral and charged states in their electronic ground states (monocation: K1, monoanion: A1, bianion: A2, trianion: A3; and tetraanion: A4) were optimized at the semiempirical AM1 and DFT levels. The polarized split-valence SVP16 basis sets and B3LYP17 functional were used for the DFT geometry

10.1021/jp906079p CCC: $40.75  2009 American Chemical Society Published on Web 10/19/2009

Study of Empty C2(3)-C82 Fullerene

J. Phys. Chem. C, Vol. 113, No. 45, 2009 19659

Figure 1. Optimized AM1 structure of C2(3)-C82 in the neutral state (a); Schlegel diagram with bond (b) and atom (c) numbering schemes.

optimizations. For all static calculations the geometry was of C2 symmetry. All quantum chemical calculations were done with the Gaussian 03 program package.18 The classical molecular dynamics simulations based on the AM1 Hamiltonian have been carried out with use of the Newton-X code.19,20 The Velocity-Verlet scheme21 was used to integrate the equation of motions. Using 0.25 fs time steps, the dynamics was equilibrated for 20 ps and then it has been running for a further 40 ps. Temperature (300 K) was controlled by an Andersen thermostat.22 Absorption spectra simulations for neutral and charged states were performed by a first sampling of the configuration space with an ensemble of thermally equilibrated nuclear geometries. The geometries were generated by the dynamics simulations mentioned above. For these samples, 80 vertical transition energies and oscillator strengths between initial and final state were calculated by the semiempirical ZINDO (Zerner’s spectroscopic parametrization for intermediate neglect of differential overlap) method.23 These quantities were used to compute the appropriate Einstein coefficient B. The Gaussian function of the energy was attributed to each point, with the height corresponding to the B coefficient

and the width equal to a phenomenological broadening constant of 0.05 eV. The sum of the evaluated Gaussian functions plotted against the transition energy gives a simulated theoretical spectrum. 3. Experimental Details The C2(3)-C82 fullerene was synthesized and separated according to the scheme described in ref 6. A PG 284 potentiostat (HEKA, Germany) was used for in situ spectroelectrochemical studies in 1,2-dichlorobenzene (o-DCB, anhydrous, 99% Aldrich) solutions with tetrabutylammonium hexafluorophosphate (TBAPF6, Fluka, dried under reduced pressure at 340 K for 24 h prior to use) as supporting electrolyte. The optical spectra were measured by the UV/vis/NIR spectrometer system TIDAS (J&M, Aalen, Germany) in o-DCB. The detailed description of experimental conditions is given in ref 6. 4. Results and Discussion As indicated above, the fullerene structure has C2 symmetry. Therefore, only one-half of the bonds and atoms were selected

19660

Sˇolc et al.

J. Phys. Chem. C, Vol. 113, No. 45, 2009

Figure 2. The optimal bond lengths calculated for the neutral state of C2(3)-C82 calculated at the AM1 and DFT theoretical levels. For bond numbering see Figure 1b.

in the numbering schemes of Figure 1. The optimal geometry of the neutral state obtained from the DFT and AM1 calculations is depicted in Figure 2. The presented bond lengths are changed in intervals from 1.37 to 1.47 Å (for DFT) and 1.36 to 1.47 Å (for AM1). Both methods indicate the smallest bond lengths for position numbers 8, 26, 29, 40, 42, 44, and 54. On the other hand, the major part of the distances is over the values 1.45 Å (for DFT) and 1.46 Å (for AM1). The electrochemical reduction or oxidation leads to structural changes of the fullerene in its charged state. The comparison of bond length differences between the neutral and charged states (K1 and from A1 to A4) is given in Figure 3. As seen, the monocationic and monoanionic state of the C82 molecule provides a minimal number of bonds where the distances are changed with respect to the neutral state. The semiempirical AM1 method compared to the DFT method reveals the higher changes in bond length numbers 35, 38, 45, 47, and 56 (up to 0.02 Å). In the case of the next negative charged states we can see that both theoretical approaches provide results that show larger changes in the bond length differences with an increase of total charge. For example, these changes are starting from outer bonds 2 to 11 for the A2 state (see Figure 1b) and for A4 are spreading to inner bonds of the molecular skeleton. The most significant bond elongation for A4 exists for bonds 4, 5, 7, 9, 47, 56, 61, and 62 (ca. 0.035 Å). On the other hand, the length for bonds 2, 6, 8, 10, 18, 24, 26, 42, 52, 58, 63, and 64 are reduced (ca. 0.03 Å). The semiempirical AM1 method has the tendency to slightly overestimate the bond length changes in comparison with DFT calculations (e.g., ca. 0.01 Å for bonds 5-9). Nevertheless, the mutual general trends are comparable. For molecules with a large number of bonds, it is appropriate to define a parameter that can characterize bond changes in one number. One of the options is the bond length alternation (BLA) parameter,24 which can help us to describe bond length changes in the whole molecule upon charging. This parameter can be calculated as M

BLA )

∑ |di - dj |/M

(1)

i)1

where dj is mean value of all bond length distances obtained from all di bond lengths and M stands for the number of bonds. With respect to this definition, a very small BLA value indicates

effective aromatic structure, e.g. BLA of benzene is 0. The data collected in Table 1 show the highest value for the neutral state. The reason for the decrease of BLA and consequently equalizing of the majority of bond lengths in charged structures could be explained as a consequence of Coulombic repulsion. The abundance of the negative electric charge or positive in the charged states leads to the increase of the distances between the atoms in the whole molecule in comparison with the distances in the neutral states (see Figure 3). For the higher charged states, the distances in pentagons and hexagons reach the uniform lengths and thus the corresponding BLA value is decreased. The positive monocharging leads to smaller changes in BLA than for negative monocharging. Both theoretical approaches reveal that the increase of the total negative charge is responsible for the monotonic decrease of BLA values. The drop for A4 is 50% for DFT and 40% for AM1 in comparison with the neutral state N. The electrochemical charging causes the changes in the distribution of partial charges on different atoms in the molecule. For all charged states including the neutral state the positive B3LYP charges are mostly located on atoms 8, 9, 12, 16, 20, 23, 25, 30, 35, 37, and 59-63 in area between 0.016e and 0.081e (see Figure 4). On the other hand, the negative charge occurs generally on atoms 2, 4, 13, 15, 17, 26, 34, 36, 64, and 65 ranging from -0.033e to -0.071e. If the total charge is increased, the negative partial charges are not regularly distributed over the major part of atoms mentioned above. For example, the mean value calculated only from the negative partial charges on 70 atoms in the molecule for the A4 state is -0.0607e. The mean value calculated from the rest of the 12 positive partial charges is 0.2504 e. The presented calculations indicate that all anionic states of the fullerene can contain atoms (e.g., 7, 8, 13, and 20) with positive partial charge (up to 0.05e). The above-mentioned nonuniform distributions of the atomic partial charges for neutral and charged states are similarly reflected in the molecular electrostatic potential. As has been shown in Figure 1S (Supporting Information), the 3-dimensional mapped electrostatic potentials of the C82 cage exhibit different characteristics depending on the total electric charge. The distribution of the unpaired electron in the molecule is connected with spin densities which are defined as the differences between the R- and β-electron densities. In Figure 5, the B3LYP Mulliken spin densities are given for selected carbon atoms of the corresponding ion radicals (see Figure 1c for the atom numbering notation). The positive part of spin density for the K1 state is primarily located on atoms 3, 11, 21, 31, 44, and 45 with values ranging from 0.056e to 0.089e. The negative part is located mainly on atoms 6, 7, 15, 16, and 25 with values in the range of -0.020e and -0.031e. The anionic state A1 has the free electron localization mainly on atoms 4, 6, 16, and 36 in the range of 0.048e to 0.070e. The spin density of the A3 state is partially similar to the spin density distribution of the cationic state K1, as indicated in Figure 5 for atoms 1, 2, 9, 24, 31, and 37-45. Moreover as in the case of the anionic state the negative part of the spin density has no meaningful values which are in the area of 0.021e. The positive part is mainly located on atoms 28, 31, 34, 44, and 45 ranging from 0.053e to 0.068e. A significant impact on the optical properties has the motion of atoms in the molecule as a part of the vibrations. Our MD simulations indicate that the shape of the fullerene ball is moreor-less periodically changed around the most probable geometry. As seen in Figure 6 the histograms of BLA values possess the Gaussian distribution for the states studied. In accordance with

Study of Empty C2(3)-C82 Fullerene

J. Phys. Chem. C, Vol. 113, No. 45, 2009 19661

Figure 3. The bond length differences between the optimal neutral and charged states of C2(3)-C82 calculated at the AM1 (solid lines) and DFT (penciled bars) theoretical levels. For bond numbering see Figure 1b.

TABLE 1: BLA Values (in Å) for the Optimal Static Geometries of the Studied Molecule in Neutral and Charged States method state

DFT

AM1

N K1 A1 A2 A3 A4

0.026 0.020 0.018 0.017 0.015 0.014

0.028 0.028 0.026 0.024 0.022 0.020

the static calculations (see Table 1), the BLA values obtained for the MD geometries are shifted to the lower values for the charged states. However, the MD simulations at room temperature show that the positions of the probability maximum are shifted about 0.005 Å to higher values. The lowest overlap (ca. 25%) between the histograms of neutral state N and charged states is observed for A3 and A4. This fact indicates that the geometry changes upon the charging or discharging cycle can have different relaxation times for the A1, A2 compared to the A3, A4 states. The relevance of the interpretation and identification of the experimental absorption peaks in the UV/vis/NIR spectra are supported by the theoretical spectra given in Figure 7. For example, the simulations for the neutral state N predict six absorption bands with increasing intensities between 0.5 and 3.0 eV. The absorption pattern of fullerenes in this visible range is mainly due to the π-π* excitations and shows a remarkable structural sensitivity. The physical origin of the individual bands might be analyzed from the static calculations of ZINDO spectroscopic properties for the optimal AM1 geometry. The HOMO to LUMO transition contributes with 69% to the lowest energy peak, which is located at 1.42 eV. The next simulated peak with a maximum at 1.86 eV originates from the HOMO-1 to LUMO transition. Finally, the transition of HOMO to LUMO+1 dominates for the third simulated peak. A good correlation between our absorption spectra simulations and the experimental spectra gives additional confirmation of the C2(3)C82 molecular structure as the major C82 isomer. The electrochemical oxidation or reduction of the fullerene under non-thin-layer conditions leads to the generation of the charged states which are mixed with neutral state in bulk.

Therefore the experimental spectrum is superimposed from the absorptions of the charged and the neutral states. Three absorption maxima at 1.36, 1.26, and 0.6 eV were detected during the oxidation of the C2(3)-C82 isomer in the first oxidation step. As shown earlier,5 they can be clearly assigned to the C82 radical cation. The occurrence of further bands in the theoretical spectrum of the K1 state can explain the diffuse shape of the experimental spectrum between 1.0 to 2.0 eV. The analysis of calculated ZINDO transitions for optimal AM1 K1 geometry reveals that the single occupied molecular orbital (SOMO) contributes significantly to the first two lowest optical transitions with relatively large oscillator strengths. The transition HOMO to SOMO (0.69 eV) contributes by more than 80% . The reduction of C2(3)-C82 in the first step causes the appearance of well-defined vis/NIR bands with maxima at 1.53, 1.39, 1.09, and 0.65 eV (Figure 7, see the second row for the A1 state). The simulated bands correspond very well with the experimental data. Despite this, a set of transitions with uniformly small intensities are predicted between 1.6 and 2.5 eV. In comparison with the K1 system, the transitions from SOMO to LUMO (0.43 eV) or SOMO to LUMO+1 (0.55 eV) dominate in the first simulated absorption peak, too. When the potential scan was extended to the second reduction step (Figure 7, second row for the A2 state), new bands with maxima at 1.47, 0.76, and 0.63 eV were observed in the experimental spectrum. Our static ZINDO calculations showed that these bands originate from the HOMO to LUMO (0.64 eV), LUMO+1 (0.84 eV), or LUMO+2 (1.28 eV) transitions. They can be assigned to the C82 dianion due to the correlation of their intensities with the charge transferred in the second reduction step. The simulated spectra for the A2 state reflected well the above-mentioned experimental bands. In addition, the next welldistinguished absorption band was obtained at 2.23 eV. Due to the presence of the neutral form in solution, this band is hidden in experimental records. For the third reduction step, an increase of the two absorption maxima at 0.94 and 0.77 eV was experimentally observed. Although the maximum at 0.77 eV overlaps strongly with the dianion band at 0.76 eV, the presence of this peak in the theoretical spectrum for the A3 state allows us to assign the experimental 0.77 eV band to the trianion (see Figure 6, bottom). The static ZINDO//AM1 calculations showed that this peak originates mostly from the SOMO to LUMO+1 transition (0.77

19662

J. Phys. Chem. C, Vol. 113, No. 45, 2009

Sˇolc et al.

Figure 4. B3LYP Mulliken charges in atomic units for the optimal neutral and charged states of C2(3)-C82. For atom numbering see Figure 1c.

Figure 5. B3LYP spin densities for the optimal neutral and charged states of C2(3)-C82. For atom numbering see Figure 1c.

Figure 7. The experimental (solid lines) and theoretical (dotted lines) absorption spectra of the neutral and charged states of C2(3)-C82 obtained from the ZINDO//AM1/MD simulations. The arrows indicate the ZINDO//AM1 transitions contributing more than 55% (H ) HOMO, S ) SOMO, L ) LUMO). Figure 6. BLA histograms for the neutral and charged states of C2(3)C82 obtained from AM1/MD simulations.

eV). Moving to even more negative potentials, including the fourth reduction step, the decrease of the most intense trianion band at 0.74 eV was experimentally observed and two new

bands of a tetraanion with maxima at 1.71 and 1.08 eV emerged (Figure 7, bottom). The theoretical spectrum for the A4 state indicates the existence of a new band at 1.14 eV, which is hidden in the mixture of various redox states of the fullerene measured. The static ZINDO//AM1 calculations showed that this band is mostly created from the HOMO to LUMO+1 (1.08 eV) and

Study of Empty C2(3)-C82 Fullerene HOMO-1 to LUMO (1.23 eV) transitions. They contribute at a significance of 55% and 58%, respectively. 5. Conclusions A systematic and comparative theoretical study based on the DFT and semiempirical calculations on the neutral and charged states of the C2(3)-C82 isomer is given. The results clarify the role of charging on the static as well as the dynamic structure of the fullerene. Minimal bond changes are reported for the monocationic and monoanionic states while the most significant are found for the tetraanionic A4 state. If the total charge at the fullerene cage is increased, the negative partial charges are regularly distributed over the major part of atoms. However, all anionic states of the C82 fullerene contain certain atoms with a positive partial charge. Molecular dynamics simulations showed that the shape of the C82 fullerene ball is more or less periodically changed around the most probable geometry. The evaluated BLA histograms reveal the similarities in geometries between neutral and charged states. The smallest mutual overlaps (ca. 25.5%) for these BLA histograms were found between the neutral state N and the charged states A3 and A4. On the basis of the collected MD geometries for charged states, the absorption spectra were simulated. As these simulations account for the natural broadening of the absorption bands due to the vibronic motion, we were able to assign more reliably the individual bands in spectroelectrochemical records where a mixture of oxidation (reduction) products exists. From the methodological point of view, this work represents the initial study for the investigation of dynamical and relaxation processes of different fullerenes with a potential for the application in optoelectronic devices, particularly in photovoltaic cells. We hope that these combined experimental and molecular dynamics techniques will help us to understand the structure-properties relations of threedimensional π-conjugation with respect to the electrochemical reduction or oxidation. Acknowledgment. Financial support by the Alexander von Humboldt Foundation (Project 3 Fokoop DEU/1063827) is gratefully acknowledged. This work was also supported by the Slovak Scientific Grant Agency (Project Nos. 1/0137/09, 1/0774/ 08, and 1/0018/09) and by the Slovak Ministry of Education (Project MVTP 2007). R.Sˇ. thanks the John von Neumann Institut for the use of the supercomputer at the FZ Ju¨lich. Supporting Information Available: The visualization of the optimal AM1 structures and 3-dimensional mapped AM1 electrostatic potentials. This material is available free of charge via the Internet at http://pubs.acs.org.

J. Phys. Chem. C, Vol. 113, No. 45, 2009 19663 References and Notes (1) Kroto, H. W.; Heath, J. R.; O’Brien, S. C.; Curl, R. F.; Smalley, R. E. Nature 1985, 318, 162. (2) Kra¨tschmer, W.; Lamb, L. D.; Fostiropoulos, K.; Huffman, D. R. Nature 1990, 347, 354. (3) Fowler, P.; Manolopoulos, D. E. An Atlas of Fullerenes; Clarendon: Oxford, UK, 1995. (4) Kikuchi, K.; Nakahara, N.; Wakabayashi, T.; Suzuki, S.; Shiromaru, H.; Miyake, Y.; Saito, K.; Ikemoto, I.; Kainosho, M.; Achiba, Y. Nature 1992, 357, 142–145. (5) Zalibera, M.; Rapta, P.; Dunsch, L. Electrochem. Commun. 2007, 9 (12), 2843–2847. (6) Zalibera, M.; Popov, A. A.; Kalbacˇ, M.; Rapta, P.; Dunsch, L. Chem.sEur. J. 2008, 14, 9960. (7) Casida, M. E.; Jamorski, C.; Casida, K. C.; Salahub, D. R. J. Chem. Phys. 1998, 108, 4439. (8) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from algorithms to applications; Academic Press: San Diego, CA, 2000. (9) Knospe, O.; Jungwirth, P. Chem. Phys. Lett. 2000, 317, 529. (10) Bredow, T.; Jug, K. Theor. Chem. Acc. 2005, 113, 1. (11) Elliotta, J. A.; Shibuta, Y. Mol. Simul. 2008, 34, 891. (12) Wang, B.-Ch.; Chen, L.; Lee, K.-J.; Cheng, Ch.-Y J. Mol. Struct. THEOCHEM 1999, 469, 127. (13) Parusel, A. B. J. J. Photochem. Photobiol. B 2000, 55, 188. (14) Scho¨nberger, H.; Schwab, Ch. H.; Hirsch, A.; Gasteiger, J. J. Mol. Model. 2004, 6, 379. (15) Dewar, M.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 2338. (16) Dunning, T. H.; Hay, P. J. In Methods of electronic structure Theory; Schaefer, H. F., III, Ed.; Plenum Press: New York, 1977; Vol. 2. (17) Stephens, P.; Devlin, P.; Chabalowski, C.; Frisch, M. J. Phys. Chem. 1994, 98, 11623. (18) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M.-C; Farkas, O.; Malick, D. K; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. GAUSSIAN 03, Revision A.1; Gaussian, Inc., Pittsburgh, PA, 2003. (19) Barbatti, M.; Granucci, G.; Ruckenbauer, M.; Persico, M.; Lischka, H. Newton-X: a package for Newtonian dynamics close to the crossing seam, version 0.13b, www.univie.ac.at/newtonx, 2006. (20) Barbatti, M.; Granucci, G.; Persico, M.; Ruckenbauer, M.; Vazdar, M.; Eckert-Maksiæ, M.; Lischka, H. J. Photochem. Photobiol. A 2007, 190, 226. (21) Verlet, L. Phys. ReV. 1968, 165, 201. (22) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (23) Zerner, M. C.; Loew, G. H.; Kirchner, R. F.; Mueller-Westerhoff, U. T. J. Am. Chem. Soc. 1980, 102, 589. (24) Jacquemin, D.; Perpete, E. A.; Chermette, H.; Ciofini, I.; Adamo, C. Chem. Phys. 1998, 79, 332.

JP906079P