Cationic Methylene–Pyrene Isomers and Isomerization Pathways

Nov 24, 2015 - News Ed., Am. Chem. .... These calculations are performed with the deMonNano code(47) using ..... calculations instead of sixth and sev...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCA

Cationic Methylene−Pyrene Isomers and Isomerization Pathways: Finite Temperature Theoretical Studies Published as part of The Journal of Physical Chemistry A virtual special issue “Spectroscopy and Dynamics of Medium-Sized Molecules and Clusters: Theory, Experiment, and Applications”. Mathias Rapacioli,†,§ Aude Simon,*,†,§ Charlotte C. M. Marshall,‡,¶,∥ Jérôme Cuny,† Damian Kokkin,‡,¶,⊥ Fernand Spiegelman,† and Christine Joblin‡,¶ †

Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France ‡ Université de Toulouse, UPS-OMP, IRAP, 31400 Toulouse, France ¶ CNRS, IRAP, 9 Avenue du Colonel Roche, BP 44346-31028 Toulouse Cedex 4, France S Supporting Information *

ABSTRACT: This paper provides spectral characterizations of the two isomers of the 1methylenepyrene cation, namely, the 1-pyrenemethylium and a pyrene-like isomer owing a tropylium cycle. Both are possible photodissociation products of the 1-methylpyrene cation and were proposed as potential contributors to the diffuse interstellar bands. In that respect, vibrational and electronic spectra are computed for the optimized structures at the density functional theory (DFT) and time-dependent (TD-)DFT levels. Finite temperature effects on these spectra are estimated from molecular dynamics simulations within the density functional-based tight-binding (DFTB) and TD-DFTB frameworks, these methods being first benchmarked against DFT and TD-DFT calculations. The computed spectra allow discrimination of the two isomers. When the temperature increases, bands are observed to redshift and merge. The isomerization mechanism is investigated with the metadynamics technique, a biased dynamics scheme allowing to probe reaction mechanisms with high energy barriers by investigating the free energy surface at various temperatures. Four pathways with similar barrier heights (3.5−4 eV) are found, showing that the interconversion process would only occur in interstellar clouds under photoactivation. The present study opens the way to simulations on larger methyl- and methylenePAHs of astrophysical interest and their experimental investigation.

1. INTRODUCTION The initial proposal of the presence of polycyclic aromatic hydrocarbons (PAHs) in the interstellar medium has motivated spectroscopic studies to identify the species that prevail in these environments.1 The astro-PAHs are the carriers of the aromatic infrared bands (AIBs), a set of emission bands observed in the 3− 14 μm range in regions where these molecules are heated by ultraviolet photons from stars.2,3 They have also been proposed as the carriers of some of the diffuse interstellar bands (DIBs), a series of absorption bands in the near-IR−visible domain observed on the extinction curve of our galaxy.4,5 Electronic spectroscopy experiments in a cold environment led Léger et al.6 to propose photoproducts of the 1-methylpyrene cation, possibly the 1-methylenepyrene cation, as carriers of four bands in the 400−850 nm region, which could match four DIBs including the strongest one at 4428 Å. Nearly 20 years later, first spectroscopic measurements were performed on gas-phase 1-methylpyrene cation and its photoproducts in the PIRENEA setup, a cold iontrap dedicated to astrochemistry.7 These experiments, which use a multiphoton-dissociation (MPD) technique, were comple© XXXX American Chemical Society

mented by calculations of the excited electronic states with the time-dependent density functional theory (TD-DFT) method.8,9 This combined approach allowed the authors to assign the 4440 and 4560 Å bands, observed in a neon matrix to the 1methylpyrene cation, while that at 4180 Å can be assigned to the 1-methylenepyrene cation, the photodissociation product resulting from the loss of a H atom from the 1-methylpyrene cation. This system (C17H+11) presents two isomeric quasidegenerate forms according to DFT calculations:7(i) the 1pyrenemethylium isomer [C16H9−CH2]+ quoted as PyrCH+2 and (ii) an isomer with a tropylium cycle quoted as PyrC+7 . Their topological formulas are reported in Figure 1. The present theoretical study is thus motivated by the fields of astrochemistry and laboratory astrophysics, more specifically: (i) The search for the carriers of the AIBs. In that respect, vibrational spectra of numerous PAHs and PAH-derived Received: September 29, 2015 Revised: November 24, 2015

A

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

tion with the BLYP functional44,45 and the 6-31G(d,p) basis set, using the aforementioned optimized geometries. This methodology was first benchmarked by Hirata and co-workers,27 who showed that it provides excitation energies for PAH radical cations such as naphthalene and anthracene consistent within 0.3 eV with the experimental data. This methodology was further applied to compute the electronic spectra of a variety of PAH+ radical cations28 and HnPAH+ molecular ions.26 DFT and TDDFT calculations are performed using the Gaussian 09 package.46 The finite temperature simulations are achieved with the energy computed via DFTB, an approximated DFT scheme that displays a low computational cost thanks to the use of parametrized integrals and a small valence basis set.35−40 In the present study, we use the second-order formulation of DFTB, or self-consistent-charge DFTB37 referred to as DFTB for sake of simplicity. These calculations are performed with the deMonNano code47 using the mio-set for Slater−Koster tables.37 As in previous works,20−25 finite temperature IR spectra are derived from MD trajectories, computing the Fourier transform of the autocorrelation function of the dipole moment.

Figure 1. Topological formulas of the two studied isomers of C17H+11: (a) PyrCH+2 and (b) PyrC+7 .

species have been computed at the DFT level in the harmonic approximation.10−17 These are absorption spectra at 0 K, whereas AIBs are emission IR spectra from photoexcited species. Assuming fast internal conversion to the electronic ground state, we may infer that the AIBs are emitted by species with high internal energy. In this context, there is a need for systematic computation of the finite temperature spectra of PAHs including the methylenepyrene cation. Previous studies have already been performed along these lines on a few specific PAHs and PAH-derived species.18−25 (ii) The search for the carriers of the DIBs. In that respect, many theoretical electronic spectra using the TD-DFT method have been determined for closed shell or radical molecular ions derived from PAHs over the past 15 years.26−34 However, experimental studies to determine electronic spectra using, for instance, the MPD technique involve finite temperature conditions. To provide theoretical background for this type of experiment, the determination of finite temperature electronic spectra is necessary. (iii) The need to determine the mechanisms/energetics of processes to ascertain the conditions (interstellar space or experimental) in which the reactions are feasible. In this regard, insights into the thermodynamics of the isomerization pathways are suitable. The goal of the present paper is thus to propose spectroscopic signatures that discriminate the isomers of 1-methylenepyrene cation in the IR and UV−visible domains, which may motivate further experimental studies. Technically, we use the DFT and TD-DFT approaches to compute vibrational and electronic spectra at the optimized ground-state geometries. Temperature effects are investigated using the density functional-based tightbinding (DFTB) and time-dependent (TD-)DFTB41 schemes coupled with molecular dynamics (MD) simulations. DFTB and TD-DFTB approaches are benchmarked against DFT and TDDFT results at T = 0 K. Finally, the free energy surfaces (FES) of the isomerization pathways for the PyrCH 2+ to PyrC 7+ interconversion in the ground-state are probed using the metadynamics technique. The computational details are reported in Section 2. Sections 3 and 4 are dedicated to present and discuss the results regarding vibrational and electronic spectra at 0 K and at finite temperature. Finally, energetics and isomerization pathway investigations are reported in Section 5.

α(ω) ∝ ω 2

∫0

+∞

dt ⟨μ(0) ·μ(t )⟩e iωt

(1)

where μ(t) is the dipole moment vector at time t, and ⟨μ(0)·μ(t)⟩ is a statistical average. The system is first thermalized at the desired temperature by means of short simulations (1 ps) in the canonical ensemble using a Nose−Hoover chain of thermostats.48−50 This dynamic allows the generation of initial conditions for MD simulations performed in the microcanonical ensemble, from which the molecular dipole moments are derived. This was done for both isomers at each studied temperature. The convergence of the band positions at 10 K was easily obtained, whereas intensity convergence required longer simulations: a fair agreement between the T = 10 K spectrum and the harmonic spectrum in terms of relative intensities was obtained for a total of 2 ns simulation (average of 20 simulations of 100 ps with a time step of 0.1 fs). To generate the temperature-dependent electronic absorption spectra, we used the TD-DFTB framework41 for computing electronic transitions. A series of 10 MD simulations on the DFTB ground state of 20 ps each with a time step of 0.2 fs were conducted in the canonical ensemble with a Nose−Hoover chain of seven thermostats with a frequency of 800 cm−1. TD-DFTB calculations (104 calculations) were then performed for snapshots extracted from these trajectories, and all the absorption spectra were summed to build a histogram. The isomerization mechanism is investigated through metadynamics simulations.51−53 These are achieved in the same ensemble and with the same thermostat characteristics as the temperature-dependent TD-DFTB simulations. The integration time step is set to 0.2 fs. In metadynamics, the MD simulation is biased to favor a specific reaction channel adding Gaussian functions to the potential energy surface (PES). In the present work, these Gaussian functions are defined by height and width of 0.027 eV and 0.1 Å, respectively. They were added every 100 MD steps, and we used the well-tempered formulation of metadynamics with a bias factor of 70.54 All metadynamics simulations were performed with the PLUMED 1.2 package55 in combination with the deMonNano code.47

2. COMPUTATIONAL DETAILS In the present study, the electronic structure is determined using two levels of theory, namely, DFT and DFTB. DFT geometry optimizations and harmonic frequency calculations are performed with the B3LYP exchange-correlation functional42 in conjunction with the 6-31G(d,p) basis set,43 while TD-DFT calculations are performed in the Franck−Condon approximaB

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

3. VIBRATIONAL SPECTRA The DFT and DFTB harmonic spectra of the two isomers are reported in Figure 2. The positions and intensities of the most

shifts for a given mode, due to external perturbations, are in fair agreement with experimental data, which demonstrates the relevance of the approach.23,56 Thus, despite these discrepancies, the DFTB and DFT spectra present overall similar main features, and both allow for the discrimination of the isomers, as discussed hereafter. In the soft-mode region (0−700 cm−1), the band corresponding to the CHdef mode (out-of-plane deformation involving CH bonds) is more intense for PyrC+7 , in particular, in terms of relative intensities. In the γCH region (∼800 cm−1), the spectrum of PyrCH+2 has additional bands with respect to PyrC+7 . These are close in energy in the DFTB spectra (see Table 1). The lowtemperature spectra differ, but, as shown below, when the temperature increases, these bands merge, making the distinction less trivial. For the in-plane bending δCH and stretching νCC mode regions (1200−1800 cm−1), the νCC modes are more intense and higher in energy for PyrCH+2 , whereas the lower-energy δCH modes are more intense for PyrC+7 . We now discuss the finite-temperature IR spectra derived from MD simulations. The spectra for PyrCH+2 and PyrC+7 were computed for 10 temperatures in the 10 to 2000 K range. A selection of spectra is reported in Figure 3. As expected, isolated bands (uncoupled modes) are redshifted when T increases, whereas bands closely lying and representing coupled/similar modes merge. In the present case, the latter are located in the δCH−νCC region (∼1200−1700 cm−1). When increasing the temperature, the spectra of the two isomers remain different, particularly in the [δCH−νCC] region. The positions of the most intense bands in the CHdef, γCH, δCH, νCC, and νCH regions as a function of temperature were computed. Those for CHdef, γCH, and νCC modes are reported in Figure 4. Those for δCH and νCH are reported in Figure S1 of the Supporting Information. For all these modes, linear behaviors are observed. The slopes of the linear fits, representing the anharmonicity factors, were determined and are reported in Table 1 (except for the νCH mode of PyrC+7 as two bands merge into one when the temperature increases). The χ values determined for the γCH modes of PyrCH+2 and PyrC+7 (present work) are found to be smaller than those for C16H10, C16H+10, and [SiC16H10]+ (1.85, 2.31, and 1.80 × 10−2 cm−1·K−1, respectively).20 We may notice that the γCH band of PyrCH+2 lying at 876 cm−1 in the harmonic approximation (Table 1) has a χ factor that is half the value of the most intense γCH band at 829 cm−1 (Table 1). Regarding the νCC mode, the χ values are different for the two isomers, that of PyrCH+2 (6.75 × 10−2 cm−1·K−1) being of the order of magnitude of the χ values for the νCC modes of −2 + C16H0/+ cm−1·K−1)20 10 and SiC16H10 (6.43/6.24 and 6.49 × 10 although slightly larger. The χ value for the νCC mode of PyrC+7 is much smaller (4.81 × 10−2 cm−1·K−1), which is more in line with that found for C24H12 (4.4 × 10−2 cm−1·K−1).20 We also estimated the χ values for the CHdef modes for the two isomers (1.32 × 10−2 and 1.41 × 10−2 cm−1·K−1), but no data is available for validation, neither for the δCH mode. Regarding the νCH mode (stretching mode at ∼3000 cm−1), the χ factors determined in the present paper are expected to be too large (no experimental data is available), as it was in the case for bare PAHs.21 Taking into account nuclear quantum effects may improve these results.

Figure 2. Harmonic spectra of PyrCH+2 (black) and PyrC+7 (red) obtained at the DFT (B3LYP) (a) and DFTB (b) levels of theory. For a given spectrum, the intensities are normalized with respect to the highest intensity band.

relevant normal modes are reported in Table 1. The positions and intensities of all normal modes are reported in Table S1 of the Supporting Information. As discussed in a previous paper,20 the absolute positions of the DFTB harmonic bands of PAH-like species differ from those of DFT by a few tens of inverse centimeters for the γCH mode (out-of-plane CH bending mode) and up to ∼200 cm−1 for the νCC (CC stretch) and νCH (CH stretch) modes. As expected, such discrepancies are observed for PyrCH+2 and PyrC+7 (cf. Table 1). However, we showed in previous studies that the position

4. ELECTRONIC SPECTRA Motivated by the longstanding issue of the assignment of the DIBs and by experimental preliminary studies,7 the electronic spectra of PyrCH+2 and PyrC+7 were also determined. TD-DFT calculations in the Franck−Condon approximation were C

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. Harmonic DFT (B3LYP) and DFTB spectra for PyrCH+2 and PyrC+7 B3LYP

DFTB

PyrCH+2 mode

σa

PyrC+7 Ia

σa

PyrCH+2 σa

Ia

Ia

PyrC+7 χ

σa

Ia

χ

204 668

7.7 21.8

1.41

CCCdef CHdef

200 732

8.9 10.0

218 703

8.3 29.3

187 699

7.6 7.5

1.32

γCH

866 887 1019

22.0 65.1 14.4

882

101.9

816 829 876

14.9 42.3 26.1

1.61 0.80

826

67.5

1.73

δCH νCC

1268 1655

104.3 395.4

1243 1534

137.7 265.9

1368 1821

48.9 272.4

1.42 6.75

∼3200

weak

∼3200

weak

∼3060

weak

7.70

71.9 202.4 167.4 weak

3.11 4.81

νCH

1362 1580 1657 ∼3060

Positions σ (in cm ) and intensities (in km·mol ) of the most relevant bands. The values of the anharmonicity factor χ (in 1 × 10−2 cm−1·K−1) derived from the finite-temperature spectra are also reported. a

−1

−1

→ S7, f = 0.22, respectively). All the aforementioned bands correspond to π → π* transitions. Table 2 also reports the 10 lowest transition energies computed at the TD-DFTB level. The general trends reported for TD-DFT spectra are well-reproduced by TD-DFTB calculations: the spectrum of PyrCH+2 is dominated by three bands between 2 and 3 eV, whereas the spectrum of PyrC+7 displays only one intense band below 3 eV and three major bands above 3.5 eV. We notice that these bands appear in the same energetic order in both calculations, except for the two bands at 3.79 and 3.97 eV for PyrC+7 , which are the eighth and ninth lines in TD-DFTB calculations instead of sixth and seventh lines at the TD-DFT level. One can observe that TD-DFTB transition energies are red-shifted with respect to TD-DFT values. Computing the electronic spectra within the TD-DFTB scheme involves the DFTB approximations, which are expected to induce additional errors to those of TD-DFT. Benchmark studies on a pool of unsaturated organic molecules are available,41 and the average deviation for adiabatic excitation energies of singlet spin states with respect to experiment was found to be 0.38 eV for TD-DFTB versus 0.36 eV for TD-DFT. We applied a uniform scaling factor of 1.075 to TD-DFTB transition energies, which led to spectra in good agreement with TD-DFT ones (see Figure 5), the maximum discrepancy being smaller than ∼0.2 eV for the most intense bands (to be compared to the 0.3 eV error usually reported between TD-DFT and experimental electronic excitation energies for PAH radical cations27). Therefore, the scaled TD-DFTB results can be considered as relevant. In the following, this scaling factor is systematically applied to transition energies. Temperature-dependent absorption spectra (TD-DFTB) are displayed in Figure 6. For both isomers, increasing the temperature results in a redshift and a broadening of all the significant bands. The shift observed between 20 and 500 K varies from 0.08 to 0.13 eV depending on the band. T = 500 K appears to be a critical temperature above which the groups of three bands observed between 2 and 3 eV for PyrCH+2 start to merge into a single broad band. The same effect is observed for the three bands between 3.5 and 4.5 eV in the PyrC+7 spectrum. At 1000 K, only the lowest transition state of PyrC+7 can be clearly identified.

Figure 3. Finite-temperature IR spectra (normalized cross-section) of PyrCH+2 (black) and PyrC+7 (red).

performed for both PyrCH+2 and PyrC+7 isomers. Twenty excited states were required for each isomer, leading to transitions from the ground-state S0 to the excited states Sx (x = 1−20). The resulting spectra in the [0−4.6 eV] range of energy are reported in Figure 5, and the 10 lowest transitions are reported in Table 2. As can be seen in Figure 5, the comparison between the electronic spectra in this energy range allows us to discriminate the two isomers. Indeed, the spectrum of the PyrCH+2 isomer is dominated by two intense bands in the [2.5−3 eV] region, whereas the PyrC+7 one presents the most intense features in the [4−4.5 eV] domain. More precisely, in the case of PyrCH+2 , three bands are found at 2.22, 2.51, and 2.92 eV (S0 → S1, S0 → S2, and S0 →S3) with oscillator strengths ( f) of 0.06, 0.12, and 0.15, respectively. Above 3 eV, only one significant band at 4.32 eV (S0 → S9) is observed. In the case of PyrC+7 , the lowest energy transition (S0 → S1, f = 0.11) is found at 2.72 eV, the next intense transition (f = 0.11) occurring at 3.64 eV (S0 → S4). The most intense bands occur at 4.08 and 4.29 eV (S0 → S6, f = 0.37 and S0 D

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 5. Electronic spectra of the two isomers obtained at the TD-DFT (no scaling) and TD-DFTB (scaling of transition energies by 1.075) levels of theory.

Table 2. TD-DFT and TD-DFTB Excitation Energies and Oscillator Strengthsa PyrCH+2 TD-DFT

TD-DFTB

ΔE

f (1 × 10−2)

ΔE

f (1 × 10−2)

2.22 2.51 2.92 3.44 3.92 3.95 4.06 4.13 4.32 4.44

6.1 12.4 15.2 1.4 0.0 3.0 0.0 0.6 6.8 2.5

2.08 2.41 2.67 2.84 2.98 3.11 3.54 3.59 3.83 4.06

4.3 16.0 13.7 0.0 0.0 1.5 0.0 2.4 0.0 4.6

PyrC+7 TD-DFT

a

TD-DFTB −2

ΔE

f (1 × 10 )

ΔE

f (1 × 10−2)

2.72 2.84 3.08 3.64 3.76 4.08 4.29 4.46 4.61 4.64

11.5 1.0 0.3 10.7 1.8 37.0 21.7 2.7 0.0 0.9

2.61 2.75 2.91 3.42 3.49 3.64 3.66 3.79 3.97 4.04

12.1 0.1 0.8 11.7 0.0 2.3 0.0 20.3 24.1 0.0

No scaling factor has been applied.

5. ENERGETICS AND ISOMERIZATION MECHANISMS In the temperature range of the simulations described in Sections 3 and 4, the simulations were not long enough to observe isomer interconversion. For simulations with temperatures higher than 2000 K, isomerization starts to be observed in competition with dissociation. To investigate at which conditions (in the interstellar medium or in the laboratory) the interconversion between the two isomers may occur, the FES for the isomerization pathways are investigated using metadynamics, a biased dynamics scheme allowing to probe a specific reaction

Figure 4. Positions of some relevant modes (see Table 1) CHdef (a), γCH (b), and νCC (c) as a function of T for PyrCH+2 (black ●) and PyrC+7 (red ◆). The data points are obtained from the temperature-dependent IR spectra; the lines result from a linear fit of the data. E

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 6. Electronic spectra of the two isomers (upper = PyrCH+2 ; lower = PyrC+7 ) as a function of temperature.

channel in a reasonable simulation length. We first checked that the DFTB approach reproduces the quasi-degeneracy of the two isomers found at the DFT level.7 The DFT and DFTB energy differences between the two isomers are reported in Table 3. At the DFT level of theory, using

PyrCH+2 taking into account ZPE effects. Thus, although PyrC+7 tends to be slightly more stable than PyrCH+2 at the DFT level, it is slightly less stable at the DFTB level. Obviously, there are various approximations in the self-consistent charge DFTB scheme (Taylor development of the energy up to second order or neglecting of three centers integrals) that could be at the origin of this difference. It is not straightforward to track the exact cause for this discrepancy, as the total energy calculation results from a self-consistent process involving the diagonalization of a matrix defined using many parameters. Revisiting the parametrization procedure of the electronic matrix elements,57 or the fitting of the repulsion term,58 might be relevant, but this is clearly beyond the scope of the present work. A difference of 8 kJ mol−1 (less than 0.1 eV) is in the range of the usually accepted DFTB errors and much smaller than the isomerization barriers (derived below). To determine the FES for the interconversion process, we used the metadynamics technique in which the MD trajectory is biased by a potential defined as a sum of Gaussian functions themselves defined according to one or several collective variables, their choice being of crucial importance to properly describe the process of interest. At the end of the simulation, the sum of all the Gaussian functions provides an approximation of the FES along the chosen collective variables. This method has the advantage of having a very low computational cost when

Table 3. DFT and DFTB Energy Differencesa (in kJ mol−1) between PyrC+7 s and PyrCH+2 ΔEel ΔEel + ΔZPE

B3LYP

BLYP

PBE0

DFTB

−4.3/−3.3 −0.5/+0.5

−5.3/−4.7 −1.5/−0.9

−3.3/−1.7 +0.6/+2.2

+4.9 +8.3

a

The results of single-point energy calculations with a larger basis set (cc-pvtz7) in the corresponding DFT/6-31G(d,p) optimized geometries are reported in italics.

the B3LYP functional, the PyrC+7 isomer was found to be 0.5 kJ mol−1 below the PyrCH+2 isomer taking into account zero-point energy (ZPE) effects (computed in the harmonic approximation). Adding ZPE corrections or increasing the size of the basis set tends to decrease the energy difference between the two isomers (see Table 3). The near-degeneracy of the two isomers was confirmed by using other functionals as seen in Table 3. At the DFTB level, the PyrC+7 isomer is found 8.3 kJ mol−1 above F

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A using DFTB for the electronic PES, therefore allowing for an exhaustive exploration of the FES. The difficulty lies on the choice of a relevant set of collective variables.

Figure 7. PyrCH+2 isomer and the corresponding labeling of the atoms.

For the specific case investigated here, starting from PyrCH+2 and using the atom labeling of Figure 7, the isomerization may proceed considering the following pathways: • a proton (H1 or H2) transfer from carbon C2 to carbon C1. A relevant coordinate for this process is dH = dC2−H1 − dC1−H1 where d stands for the distance. We noticed that the choice of proton H1 or H2 does not matter here, as H1 and H2 are regularly interchanged during the dynamics. • the insertion of carbon C2 either between carbons C1 and C3 (hereafter called the α mechanism) or between carbons C1 and C4 (β mechanism). A natural choice for the coordinate of the α mechanism is the variable dCα = dC1−C3 − dC2−C3 and dCβ = dC1−C4 − dC2−C4 for the β mechanism. To avoid either dissociations or isomerizations that are not of interest in the present study, we applied constraints to the system, consisting of potential energy walls preventing migration of hydrogen atoms to other carbon atoms or loss of cyclic patterns. Figure 8 reports the FES obtained at 300 K for both the α and β mechanisms. A and C basins correspond to the PyrCH+2 and PyrC7+ isomers, respectively. Considering first the α-type isomerization (A → C, upper map of Figure 8), one clearly sees that two distinct paths are possible. The first one corresponds to a proton transfer (A→ B) followed by a carbon insertion (B→ C), whereas the second one corresponds to a carbon insertion (A→ D) followed by a proton transfer (D→ C). Interestingly, it seems very unfavorable that the two processes occur simultaneously as a high free-energy barrier raises at the map origin. The same conclusions hold for the β-type mechanism. The identification of these four possible distinct isomerization paths yields the isomerization scheme of Figure 9, displaying seven characteristic paths (notice that the same path A→ B is present in one α and one β isomerization mechanism). For each path, we performed metadynamics simulations at various temperatures, constraining the system to explore only this specific path. Results are presented in Figure 10. It is interesting to note that temperature has a very minor effect on the free-energy barriers involved in the isomerization. This can be understood because of the moderate temperature range exploration and because minor entropic effects are expected in the case of a single atom (hydrogen) or small functional groups (CH2) migrations and ring openings. The determined isomerization pathways do not involve soft motions of large functional groups. For the four possible paths leading to isomerization

Figure 8. Free energy surfaces corresponding to the isomerization of type α (upper) and β (lower).

(P1+P2; P1+P5; P3+P4; P6+P7), the free energy barrier is ∼3.5−4 eV. This is consistent with the absence of isomerization in the temperature range of the MD simulations presented in Sections 3 and 4. Conditions for isomerization are, however, likely to be fulfilled in the MPD experiments performed by Kokkin and co-workers.7 In interstellar space, this isomerization process would occur in regions irradiated by UV photons. Given the energy range, excited states are expected to play a role in the isomerization mechanism. Although this is not the scope of the present paper, it would be worth investigating this mechanism. It is interesting to compare the barrier heights found here for the isomerization process between PyrCH+2 and PyrC+7 with those already reported at the DFT level for the isomerization from benzylium to tropylium ions.59 Indeed, although the systems are different, they display some chemical similarities suggesting that their isomerization pathways could be similar. In the benzylium case, the α and β mechanisms are equivalent due to symmetry. Bullins et al.59 investigated the process that would correspond to the P3/P6+P4/P7 paths in our system. The equivalent of isomers D/E is found to be 2.4 eV above the equivalent of A and 2.8 eV above the equivalent of C. This can be compared with values of ∼2.7(D)/2.95(E) and 2.95(D)/ 3.05(E) eV extracted from the low-temperature curves of Figure 10. We observe the same trends for the PES of the benzylium as for the FES of the 1-pyrenemethylium, despite the differences in both their chemical nature and the level of theory. However, one should note that, in the present study, an additional competitive mechanism has been identified for this isomerization (P1+P2/ P5), and it would be worth investigating the equivalent path for the benzylium/tropylium isomerization.

6. CONCLUSION The present paper provides vibrational (mid-IR) and electronic (UV−visible) characterizations of the two isomer forms of the 1methylenepyrene cation (or dehydrogenated 1-methylpyrene cation), a molecular ion of potential astrophysical interest. By G

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 9. Reaction paths for the isomerization. Paths for α-type and β-type isomerizations are represented with plain and dashed lines, respectively.

Figure 10. Free energy curves along some characteristic reaction paths (see Figure 9) computed at different temperatures.

isomers. Finite temperature spectra are reported, which should be of use both in the interpretation of laboratory MPD experiments and in the analysis of astronomical spectra (the

performing static DFT/TD-DFT and dynamic DFTB/TDDFTB dynamic calculations, we show that both vibrational and electronic spectra allows for the discrimination of the two H

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

(4) Herbig, G. H. The Diffuse Interstellar Bands. Annu. Rev. Astron. Astrophys. 1995, 33, 19−73. (5) Salama, F.; Galazutdinov, G. A.; Krełowski, J.; Biennier, L.; Beletsky, Y.; Song, I.-O. Polycyclic Aromatic Hydrocarbons and the Diffuse Interstellar Bands: A Survey. Astrophys. J. 2011, 728, 154. (6) Lég er, A.; D’Hendecourt, L.; Defourneau, D. Proposed Identification for the (Common) Carrier of the 4430 Å, and 7565 Å, DIBs. Astron. Astrophys. 1995, 293, L53−L56. (7) Kokkin, D.; Simon, A.; Marshall, C.; Bonnamy, A.; Joblin, C. A Novel Approach to the Detection and Characterization of PAH Cations and PAH-Photoproducts. In Proceedings of the International Astronomical Union, IAU Symposium No. 297, 2013; Cami, J., Cox, N. L. J., Eds.; IAU publisher, Cambridge University Press, 2014; Vol. 297, pp 286− 290.10.1017/S1743921313016001 (8) Runge, E.; Gross, E. K. U. Density Functional Theory for TimeDependent Systems. Phys. Rev. Lett. 1984, 52, 997−1000. (9) Petersilka, M.; Gossmann, U. J.; Gross, E. K. U. Excitation Energies from Time-Dependent Density-Functional Theory. Phys. Rev. Lett. 1996, 76, 1212−1215. (10) Bauschlicher, C. W., Jr.; Peeters, E.; Allamandola, L. J. The Infrared Spectra of Very Large Irregular Polycyclic Aromatic Hydrocarbons (PAHs): Observational Probes of Astronomical PAH Geometry, Size, and Charge. Astrophys. J. 2009, 697, 311−327. (11) Ricca, A.; Bauschlicher, C. W., Jr.; Mattioda, A. L.; Boersma, C.; Allamandola, L. J. The Far-Infrared Spectroscopy of Very Large Neutral Polycyclic Aromatic Hydrocarbons. Astrophys. J. 2010, 709, 42−52. (12) Ricca, A.; Bauschlicher, C. W., Jr.; Boersma, C.; Tielens, A. G. G. M.; Allamandola, L. J. The infrared Spectroscopy of Compact Polycyclic Aromatic Hydrocarbons Containing up to 384 Carbons. Astrophys. J. 2012, 754, 75. (13) Bauschlicher, C. W., Jr.; Ricca, A. The Infrared Spectra of Polycyclic Aromatic Hydrocarbons with Some or All Hydrogen Atoms Removed. Astrophys. J. 2013, 776, 102. (14) Mackie, C. J.; Peeters, E.; Bauschlicher, C. W., Jr.; Cami, J. Characterizing the Infrared Spectra of Small, Neutral, Fully Dehydrogenated Polycyclic Aromatic Hydrocarbons. Astrophys. J. 2015, 799, 131. (15) Joalland, B.; Simon, A.; Marsden, C. J.; Joblin, C. Signature of [SiPAH]+ π-Complexes in the Interstellar Medium. Astron. Astrophys. 2009, 494, 969−976. (16) Simon, A.; Joblin, C. Thermochemistry and Infrared Spectroscopy of Neutral and Cationic Iron-Polycyclic Aromatic Hydrocarbon Complexes of Astrophysical Interest: Fundamental Density Functional Theory Studies. J. Phys. Chem. A 2007, 111, 9745−9755. (17) Simon, A.; Joblin, C. The Computed Infrared Spectra of a Variety of [FePAH]+ Complexes: Mid- and Far-Infrared Features. Astrophys. J. 2010, 712, 69−77. (18) Van Oanh, N. T.; Parneix, P.; Bréchignac, P. Vibrational Dynamics of the Neutral Naphthalene Molecule from a Tight-Binding Approach. J. Phys. Chem. A 2002, 106, 10144−10151. (19) Falvo, C.; Friha, H.; Pino, T.; Dhaouadi, Z.; Parneix, P.; Calvo, F.; Brechignac, P. Effects of Hydrogen Dissociation on the Infrared Emission Spectra of Naphthalene: Theoretical Modeling. Phys. Chem. Chem. Phys. 2013, 15, 10241−10250. (20) Joalland, B.; Rapacioli, M.; Simon, A.; Joblin, C.; Marsden, C. J.; Spiegelman, F. Molecular Dynamics Simulations of Anharmonic Infrared Spectra of [SiPAH]+ π-Complexes. J. Phys. Chem. A 2010, 114, 5846−5854. (21) Simon, A.; Rapacioli, M.; Lanza, M.; Joalland, B.; Spiegelman, F. Molecular Dynamics Simulations on [FePAH]+ π-complexes of Astrophysical Interest: Anharmonic Infrared Spectroscopy. Phys. Chem. Chem. Phys. 2011, 13, 3359−3374. (22) Rapacioli, M.; Simon, A.; Dontot, L.; Spiegelman, F. Extensions of DFTB to Investigate Molecular Complexes and Clusters. Phys. Status Solidi B 2012, 249, 245−258. (23) Simon, A.; Rapacioli, M.; Mascetti, J.; Spiegelman, F. Vibrational Spectroscopy and Molecular Dynamics of Water Monomers and Dimers Adsorbed on Polycyclic Aromatic Hydrocarbons. Phys. Chem. Chem. Phys. 2012, 14, 6771−6786.

AIBs). The FES for the isomerization mechanism are determined using the metadynamics approach. Four pathways with similar barrier heights in the FES (3.5−4 eV) are identified. The heights of the barriers preclude from direct isomer interconversion on the ground state at low temperatures. In the interstellar medium, the process would require photoactivation with efficient energy funneling into the ground state with vibrational excitation. Since the energy barrier is comparable to the lowest excitation energies, electronic excited states are also likely to be involved in the isomerization process. Investigating their precise role, as achieved for other PAHs,60,61 will be of interest, but this is beyond the scope of this paper. The present work also opens the way to similar studies on larger methyl- and methylene-PAHs of astrophysical interest and provides motivation for further experimental and theoretical investigations on these systems.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b09494. Tabulated harmonic B3LYP and DFTB vibrational frequencies of for PyrCH+2 and PyrC+7 . Positions of the δCH and νCH vibrational modes plotted as a function of T for PyrCH+2 and PyrC+7 . (PDF)



AUTHOR INFORMATION

Corresponding Author

*Phone: +33561556033. E-mail: [email protected]. Present Addresses ∥

School of Chemistry, The University of Nottingham, University Park, Nottingham, NG7 2RD, U.K. ⊥ Department of Chemistry and Biochemistry, Arizona State University, Tempe, Arizona 85287, USA. Author Contributions §

Contributed equally to the work.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors acknowledge the computing facility CALMIP for allocation of computer resources and CNRS and Region MidiPyrénées. This study has been supported by the French ANR funding with reference ANR-10-BLAN-0501-GASPARIM “Gasphase PAH research for the interstellar medium” and ANR 13BS08-0005 PARCS “Polycyclic Aromatic Hydrocarbons Reactivity in Cryogenic Solids”, by the GDR Groupement de Recherche 3533 EMIE Edif ices Moléculaires Isolés et Environnés and partially supported through the NEXT grant No. ANR-10-LABX-0037 in the framework of the “Programme des Investissements d’Avenir”. C.C.M.M. thanks The Univ. of Nottingham and The British Council for financial support.



REFERENCES

(1) PAHs and the Universe: A Symposium to Celebrate the 25th Anniversary of the PAH Hypothesis; EAS Publications Series. Joblin, C., Tielens, A. G. G. M., Eds.; EAS Publications, 2011; Vol. 46. (2) Léger, A.; Puget, J. L. Identification of the Unidentified IR Emission Features of Interstellar Dust. Astron. Astrophys. 1984, 137, L5−L8. (3) Allamandola, L. J.; Tielens, A. G. G. M.; Barker, J. R. Polycyclic Aromatic-Hydrocarbons and the Unidentified Infrared-Emission Bands - Auto Exhaust Along the Milky-Way. Astrophys. J. 1985, 290, L25−L28. I

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (24) Simon, A.; Spiegelman, F. Water Clusters Adsorbed on Polycyclic Aromatic Hydrocarbons: Energetics and Conformational Dynamics. J. Chem. Phys. 2013, 138, 194309. (25) Simon, A.; Spiegelman, F. Conformational Dynamics and FiniteTemperature Infrared Spectra of the Water Octamer Adsorbed on Coronene. Comput. Theor. Chem. 2013, 1021, 54−61. (26) Hammonds, M.; Pathak, A.; Sarre, P. J. TD-DFT Calculations of Electronic Spectra of Hydrogenated Protonated Polycyclic Aromatic Hydrocarbon (PAH) Molecules: Implications for the Origin of the Diffuse Interstellar Bands? Phys. Chem. Chem. Phys. 2009, 11, 4458− 4464. (27) Hirata, S.; Lee, T. J.; Head-Gordon, M. Time-Dependent Density Functional Study on the Electronic Excitation Energies of Polycyclic Aromatic Hydrocarbon Radical Cations of Naphthalene, Anthracene, Pyrene, and Perylene. J. Chem. Phys. 1999, 111, 8904−8912. (28) Hirata, S.; Head-Gordon, M.; Szczepanski, J.; Vala, M. TimeDependent Density Functional Study of the Electronic Excited States of Polycyclic Aromatic Hydrocarbon Radical Ions. J. Phys. Chem. A 2003, 107, 4940−4951. (29) Malloci, G.; Joblin, C.; Mulas, G. Theoretical Evaluation of PAH Dication Properties. Astron. Astrophys. 2007, 462, 627−635. (30) Malloci, G.; Mulas, G.; Cecchi-Pestellini, C.; Joblin, C. Dehydrogenated Polycyclic Aromatic Hydrocarbons and UV Bump. Astron. Astrophys. 2008, 489, 1183−1187. (31) Malloci, G.; Mulas, G.; Cappellini, G.; Joblin, C. Time-Dependent Density Functional Study of the Electronic Spectra of Oligoacenes in the Charge States −1, 0, + 1, and + 2. Chem. Phys. 2007, 340, 43−58. (32) Malloci, G.; Cappellini, G.; Mulas, G.; Mattoni, A. Electronic and Optical Properties of Families of Polycyclic Aromatic Hydrocarbons: A Systematic (Time-Dependent) Density Functional Theory Study. Chem. Phys. 2011, 384, 19−27. (33) Useli-Bacchitta, F.; Bonnamy, A.; Mulas, G.; Malloci, G.; Toublanc, D.; Joblin, C. Visible Photodissociation Spectroscopy of PAH Cations and Derivatives in the PIRENEA Experiment. Chem. Phys. 2010, 371, 16−23. (34) Weisman, J. L.; Lee, T. J.; Salama, F.; Head-Gordon, M. TimeDependent Density Functional Theory Calculations of Large Compact Polycyclic Aromatic Hydrocarbon Cations: Implications for the Diffuse Interstellar Bands. Astrophys. J. 2003, 587, 256−261. (35) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Construction of Tight-Binding-like Potentials on the Basis of Density Functional Theory - Application to Carbon. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 12947−12957. (36) Seifert, G.; Porezag, D.; Frauenheim, T. Calculations of Molecules, Clusters, and Solids with a Simplified LCAO-DFT-LDA Scheme. Int. J. Quantum Chem. 1996, 58, 185−192. (37) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-binding Method for Simulations of Complex Material Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260− 7268. (38) Frauenheim, T.; Seifert, G.; Elsterner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. A Self-Consistent Charge DensityFunctional Based Tight-Binding Method for Predictive Materials Simulations in Physics, Chemistry and Biology. Phys. Status Solidi B 2000, 217, 41−62. (39) Frauenheim, T.; Seifert, G.; Elstner, M.; Niehaus, T.; Köhler, C.; Amkreutz, M.; Sternberg, M.; Hajnal, Z.; Carlo, A. D.; Suhai, S. Atomistic Simulations of Complex Materials: Ground-State and Excited-State Properties. J. Phys.: Condens. Matter 2002, 14, 3015. (40) Oliveira, A.; Seifert, G.; Heine, T.; Duarte, H. Density Functional based Tight Binding: an Approximate DFT Method. J. Braz. Chem. Soc. 2009, 20, 1193−1205. (41) Niehaus, T. A.; Suhai, S.; Della Sala, F.; Lugli, P.; Elstner, M.; Seifert, G.; Frauenheim, T. Tight-Binding Approach to Time-Dependent Density-Functional Response Theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 63, 085108. (42) Becke, A. D. Density-Functional Thermochemistry 0.3. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652.

(43) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-Consistent MolecularOrbital Methods 0.25. Supplementary Functions for Gaussian-Basis Sets. J. Chem. Phys. 1984, 80, 3265−3269. (44) Becke, A. D. Density-Functional Exchange-Energy Approximation With Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (45) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (46) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian09, Revision D.01; Gaussian Inc: Wallingford, CT, 2009. (47) Heine, T.; Rapacioli, M.; Patchkovskii, S.; Frenzel, J.; Koster, A.; Calaminici, P.; Duarte, H. A.; Escalante, S.; Flores-Moreno, R.; Goursot, A. et al. deMonNano; http://demon-nano.ups-tlse.fr/, 2009. (48) Nosé, S. A. Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511. (49) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (50) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé-Hoover Chains: The Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 1992, 97, 2635−2643. (51) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (52) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302. (53) Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the Accuracy of Metadynamics. J. Phys. Chem. B 2005, 109, 6714−6721. (54) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603. (55) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; et al. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (56) Simon, A.; Iftner, C.; Mascetti, J.; Spiegelman, F. Water Clusters in an Argon Matrix: Infrared Spectra from Molecular Dynamics Simulations with a Self-Consistent Charge Density Functional-Based Tight Binding/Force-Field Potential. J. Phys. Chem. A 2015, 119, 2449− 2467. (57) Oliveira, A. F.; Philipsen, P.; Heine, T. DFTB Parameters for the Periodic Table, Part 2: Energies and Energy Gradients from Hydrogen to Calcium. J. Chem. Theory Comput. 2015, 11, 5209−5218. (58) Wahiduzzaman, M.; Oliveira, A. F.; Philipsen, P.; Zhechkov, L.; van Lenthe, E.; Witek, H. A.; Heine, T. DFTB Parameters for the Periodic Table: Part 1, Electronic Structure. J. Chem. Theory Comput. 2013, 9, 4006−4017. (59) Bullins, K. W.; Huang, T. T. S.; Kirkby, S. J. Theoretical Investigation of the Formation of the Tropylium Ion From the Toluene Radical Cation. Int. J. Quantum Chem. 2009, 109, 1322−1327. (60) Tokmachev, A. M.; Boggio-Pasqua, M.; Bearpark, M. J.; Robb, M. A. Photostability via Sloped Conical Intersections: A Computational Study of the Pyrene Radical Cation. J. Phys. Chem. A 2008, 112, 10881− 10886. (61) Mendive-Tapia, D.; Perrier, A.; Bearpark, M. J.; Robb, M. A.; Lasorne, B.; Jacquemin, D. New Insights into the By-Product Fatigue Mechanism of the Photo-Induced Ring-Opening in Diarylethenes. Phys. Chem. Chem. Phys. 2014, 16, 18463−18471.

J

DOI: 10.1021/acs.jpca.5b09494 J. Phys. Chem. A XXXX, XXX, XXX−XXX