Computational study of the adamantane cation: simulations of spectroscopy, fragmentation dynamics, and internal conversion

Diamondoids, of which adamantane (C10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{10}$$\end{document}H16\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}$$\end{document}) is the simplest representative, constitute an intriguing class of carbon based nanomaterials with interesting chemical, mechanical and opto-electronic properties. While neutral diamondoids have been extensively studied for decades, their cationic counterparts were a subject of recent experimental investigations motivated by their potential role in astrochemistry. Here, we perform a computational study of the adamantane cation (C10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{10}$$\end{document}H16+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}^+$$\end{document}) complementing the recent experimental findings. Specifically, we extend earlier theoretical work on vibrationally resolved electronic spectroscopy by accounting for the higher lying electronically excited states of the cation in the absorption and photoelectron spectra. We also perform adiabatic and nonadiabatic (surface hopping) molecular dynamics simulations to study (fast) fragmentation processes and electronic relaxation of C10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{10}$$\end{document}H16+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}^+$$\end{document}. Our simulations reveal that after excitation with near-infrared–ultraviolet photons, the adamantane cation undergoes an ultrafast internal conversion to the ground (doublet) state (on a time scale of 10–100 fs depending on initial excitation energy) which can be followed by a fast fragmentation, predominantly H loss. The yield of the ultrafast hydrogen dissociation as a function of excitation energy is reported.


Introduction
There are different types of carbon adaptation present in nature and interstellar chemistry [1][2][3]. One interesting class of carbon compounds are the diamondoids having diamond-like, hydrogen saturated molecular structures [4]. These molecules have potential various applications in nanotechnology [5,6], biochemistry [7] as well as medicinal chemistry and pharmaceuticals industry [8][9][10][11]. They have a good chemical and thermal stability which has been a major factor of their studies in the scientific community [4,12]. The small diamondoids are naturally found in natural gases and crude oil with great abundances [13]. Nanodiamonds (large diamondoinds consisting of a few 100 carbon atoms) constitute carbonaceous, primitive chondritic meteorites [2]. And with a relatively recent detection of buckminsterfullerene (C 60 ) in the interstellar medium (ISM) [14][15][16], there have been speculations of presence of more large size carbon derivates in the ISM [17].
In addition, the EPD study by Crandall et al. [39] provides information on photoinduced fragmentation of Ada + . In particular, it was found that the H loss channel dominates at low excitation energies (with appearance energy of 1.38 eV with respect to the cationic ground state). Moreover, Candian et al. investigated the dissociation of Ada + with imaging Photoelectron Photoion Coincidence (iPEPICO) spectroscopy by exciting Ada with vacuum ultraviolet (VUV) photons [45]. The appearance energy for the H loss channel was reported to be 10.50 eV with respect to the neutral ground state (or 1.21 eV with respect to the cationic ground state) [45]. Furthermore, fragmentation of Ada + has been examined by electron ionization mass spectrometry (EIMS) [3,46]. We also note that besides fragmentation, adamantane cage-opening for the related compound, amantadine cation, was recently investigated [47,48].
In this computational study, we pursue several directions in studying Ada + . First of all, we extend our previous spectroscopy modeling [30,44] by accounting for the higher lying electronically excited states of the cation in the vibrationally resolved absorption and photoelectron spectra (albeit using a simple vertical gradient approximation as described below). Secondly, we perform molecular dynamics simulations for Ada + relying on the ground doublet state only and providing certain excess energy in order to reveal fast fragmentation processes. Thirdly, we model ultrafast nonadiabatic dynamics (internal conversion) of Ada + using the surface hopping (SH) approach.

Electronic structure methods
The adamantane neutral (Ada) and cation (Ada + ) were optimized in their ground electronic states ( S 0 and D 0 , respectively) using a series of semiempirical methods (AM1 [49], PM3 [50], PM7 [51], OM2 [52], and OM3 [53]) as well as DFT at the B3LYP [54, 55]/6-31G* [56,57] level. The spin-restricted (R) formalism was used for the neutral, whereas the cation geometry was optimized using the spinunrestricted (U) formalism. In addition, the cation geometry was optimized at the semiempirical configuration interaction singles and doubles (CISD) level, which was used for surface hopping simulations (see below). In particular, the molecular orbitals for the CISD calculation were obtained from a restricted open-shell (RO) OM3 calculation, and the CISD calculation was performed in the restricted active space of 9 occupied and 10 virtual orbitals. We term this method ROOM3/CISD(9 × 10 ). For all optimized structures, the normal mode analysis was executed to confirm the minimum nature of the found stationary points (i.e., the absence of imaginary frequencies).

Electronic spectra
The vertical excitation energies and the corresponding oscillator strengths of the D 0 → D i ( i > 0 ) electronic transitions of the cation were computed using linear response time-dependent (TD) DFT, at the TD-UB3LYP/6-31G* level, as well as with the ROOM3/CISD(9 × 10 ) method. These calculations were done at the ground-state ( D 0 ) geometries of the cation optimized with UB3LYP/6-31G* and ROOM3/CISD(9 × 10 ). In addition, the D 0 → D i vertical transitions were calculated at the optimized geometry of the neutral ( S 0 minimum), at the TD-DFT level, to estimate transition energies of a vertical photoelectron spectrum (by adding the vertical, S 0 → D 0 ionization potential to the D 0 → D i vertical excitation energies).
Further, in order to compute a vibrationally resolved absorption specturm of the cation as well as a vibrationally resolved photoelectron spectrum, we tried to optimize excited states of the cation, but we were unable to locate excited-state minima. Therefore, we resorted to the vertical gradient (VG) approximation [58], which constructs the final state ( D i ) potential energy surface (PES) by shifting the harmonic PES of the initial state ( D 0 in the case of the cation absorption and S 0 for the photoelectron spectrum) based on the gradient of the final state at the minimum geometry of the initial state. For the S 0 → D 0 transition, we used in addition the adiabatic Hessian (AH) method, which requires normal modes of the initial and final states, which, in turn, we were able to obtain due to successful optimizations for the both, S 0 and D 0 , states. The photoelectron spectra were computed using a time-independent (TI) formalism (based on the calculations of the Franck-Condon (FC) factors), and the absorption spectra using a time-dependent (TD) formalism (based on the calculations of the autocorrelation functions), since the TI formalism resulted in low spectrum progression for some of the considered transitions [58]. All vibrationally resolved spectra were computed in the FC approximation. The default Gaussian 16 parameters were used. In particular, a Gaussian broadening with half-width at half-maximum of 135 cm −1 was applied.

Born-Oppenheimer and surface hopping molecular dynamics
Further, we performed ab initio molecular dynamics simulations for the adamantane cation with the focus on possible fragmentation. To do so, we considered three scenarios, shown in Fig. 1 and labeled there (a), (b), and (c), respectively. Scenarios (a) and (b) correspond to EPD and VUV (EIMS) techniques, respectively, both under the approximation that omits explicit treatment of cationic excited states, and, thus, does not explicitly describe the process of internal conversion (IC). In (a), it is assumed that the cation, which is originally in the D 0 state, is excited with an UV (or visible, or NIR) laser and the photon energy E h is transferred to the nuclear kinetic energy E kin , and that the cation remains in the ground state ( D 0 ). In (b), it is assumed that the neutral, which is originally in the S 0 state, is excited either with VUV photons or an electron beam, and a part of excitation energy E h is needed for ionization (ionization potential) whereas the rest goes to the kinetic energy of nuclei, and the cation is in its ground state ( D 0 ). The kinetic energy of an outgoing electron is assumed to be zero (as in threshold photoelectron photoion coincidence (TPEPICO) spectroscopy). Then Born-Oppenheimer molecular dynamics (BOMD) simulations are performed in the cationic ground state. Because of energy-conservation problems found in test calculations employing unrestricted methods, Fermi smearing is used (which corresponds to some extent to an electronic excitation). In the Fermi smearing approach, the molecular orbital occupations are allowed to be fractional and are calculated using the Fermi-Dirac distribution, with an electronic temperature being a parameter [59]. The default (in the MNDO program [60]) value of 20,000 K was used in our simulations. We note that methods (a) and (b) are similar to Grimme's QCEIMS (quantum chemistry EIMS) [61] approach. In our work, the UOM3 method is used in dynamics simulations, since it was found to yield the adamantane ionization potential in good agreement with experiment (see below). The trajectories were propagated for 10 ps with a time step of 0.5 fs using velocity Verlet algorithm. The overall translation and rotation of the molecule were eliminated to focus only on the internal dynamics. The trajectories were propagated at constant total energy.
In scenario (c), which corresponds to EPD (similarly to (a)), the cationic excited states are taken into account explicitly, and the nonadiabatic dynamics after photoexcitation are simulated. Thus, IC is directly accounted for in method (c). Specifically, Tully's surface hopping (SH) approach is used here to perform nonadiabatic dynamics simulations [62]. The electronic structure of the cation is described using the semiempirical ROOM3/CISD(9 × 10 ) method as introduced above. This method yields spin-pure doublets and does not suffer from a spin-contamination problem. Moreover, the double excitations are taken into account in addition to singles in order to avoid potential problems in the description of D 0 ∕D 1 conical intersections. The 9 × 10 active space was chosen based on test calculations of excitation energies  with different active spaces. The SH trajectories were propagated for 1 ps with a nuclear time step of 0.1 fs. The energy-based decoherence correction of Granucci and Persico was applied [63].
The initial conditions for the dynamics simulations were prepared in the following way. For method (a), we performed ground-state ( D 0 ) molecular dynamics at a constant temperature of T = 20 K -the temperature given as an upper limit in the experiment by Crandall et al. [39]. For method (b), we used the same approach, but the dynamics were modeled in the S 0 state and with T = 500 K, a common temperature in mass spectrometry experiments [61]. The temperature was controlled by a simple velocity-rescaling algorithm, i.e., multiplying velocities by √ Here, T ref is the reference (desired) temperature and T curr is the current (instantaneous) temperature. The initial velocities were randomly sampled from a uniform distribution and scaled to yield the reference temperature. The trajectories were propagated for 100 ps with a time step of 0.5 fs. Then, 1000 snapshots were selected from these ground-state trajectories, starting after 50 ps and sampling every 50 fs. Only selected ground-state geometries (and not ground-state velocities) were used for the subsequent photodynamics simulations. The initial velocities (for the photodynamics simulations) were instead generated randomly and scaled to yield the desired kinetic energy E kin (see Fig. 1). For method (c), the initial conditions (geometries and velocities) were prepared by sampling a Wigner function (at 0 K) of the harmonic ROOM3/ CISD(9 × 10 ) ground state of the cation. In addition we tested these initial conditions for method (a), rescaling the Wigner velocities to yield the desired kinetic energy. Further, the vertical absorption spectra were computed for the sampled geometries using the ROOM3/CISD(9 × 10 ) method. After that, the SH simulations were performed. Four sets of SH simulations, differing by initial excitation energy, were considered: (i) excitation to the D 1 ∕D 2 states, (ii) excitation to the D 3 ∕D 4 states, (iii) excitation to an energy window (2.3,2.9) eV, and (iv) excitation to the D 8 state. For cases (i)-(iii), the initial electronic state of an individual trajectory was determined by comparing oscillator strengths of respective transitions and selecting the state corresponding to the transition with the largest oscillator strength. Each of four batches of the SH simulations included about 1000 trajectories.

Programs
The calculations at the (TD-)DFT level as well as those with PM7 were done with Gaussian 16 [64]. The calculations involving the OMx (x = 2 or 3) methods (including the SH simulations) as well as AM1 and PM3 were performed with the MNDO program [60]. The adiabatic dynamics simulations were performed with our in-house code coupled with MNDO. The calculations of the vibrationally-resolved spectra were done using Gaussian 16.

Geometries and ionization potentials
We first optimized the adamantane neutral and cation using semiempirical methods (OM2, OM3, AM1, PM3, and PM7) and DFT at the (U)B3LYP/6-31G* level. Selected geometric parameters (see Fig. 2) are shown in Table S1. The CC bond lengths of the neutral R 1 , R 2 , R 3 are equal to each other and are ∼1.53 Å at the AM1, PM3, OM2, and B3LYP levels and ∼1.54 Å at the OM3 and PM7 levels. The CH bond lengths r 1 , r 2 , r 3 , r 4 , r 5 vary between ∼1.10 Å (B3LYP) and ∼1.14 Å (OM3). We also note that AM1, PM3, and B3LYP predict r 1 , r 2 , r 3 , r 4 , r 5 to be virtually equal to each other, whereas OM2, OM3, and PM7 show deviations of up to 0.02 Å. For the cation, all the methods predict the R 2 elongation in comparison to the neutral, which is caused by the Jahn-Teller effect [65]. Specifically, R 2 ≈ 1.57 Å for AM1 and OM2, R 2 ≈ 1.58 Å for PM3 and OM3, R 2 ≈ 1.60 Å for PM7, and R 2 ≈ 1.61 Å for UB3LYP. The CH r 1 bond is also elongated in comparison to r 1 of the neutral and also in comparison to all other CH r i bonds of the cation. Overall, the geometries obtained with the used semiempirical methods are mostly consistent with current and previous DFT calculations [30,65] as well as experimental results [66,67]. Further, we calculated vertical (Vert IP), adiabatic (Adia IP) and 0-0 (also known as Fig. 2 The molecular structure of adamantane with labeling of selected bonds and angles. The picture is adapted from [30]. Copyright 2018, American Institute of Physics zero-point-energy(ZPE)-corrected adiabatic, ZPE+) ionization potentials of adamantane from the ground state ( S 0 ) of the neutral to the ground state ( D 0 ) of the cation (see Table 1). The evaluated experimental IP value from the NIST database is 9.25 ± 0.04 eV [68]. Comparing the 0-0 IPs to this value, we find that PM7 performs best ( IP 0-0 = 9.23 eV), followed by OM3 ( IP 0-0 = 9.15 eV). The B3LYP value is too small ( IP 0-0 = 8.71 eV) in agreement with previous reports [23,30]. However, we select the OM3 method for adiabatic ( D 0 ) dynamics simulations since PM7 is not implemented in the MNDO software (which, in turn, is used for our inhouse adiabatic dynamics code) and, moreover, the difference between OM3 and PM7 0-0 IPs is only 0.08 eV.

Excited states of the cation
We calculated the vertical excitation energies and oscillator strengths of the adamantane cation using TD-UB3LYP/6-31G* and ROOM3/CISD(9 × 10 ) methods. Each method was applied to both UB3LYP/6-31G* and ROOM3/CISD(9 × 10 ) optimized geometries of the cation ( D 0 minima). The excitation energies and oscillator strengths are collected in Table 2.
At the DFT geometry, the lowest two states D 1 and D 2 , which are degenerate, are located at ∼1.07 eV. At the OM3/CISD(9 × 10 ) geometry, they are found at ∼0.9 eV. The oscillator strengths of the transitions to these states are very small; at that, the semiempirical oscillator strengths are ∼4.3 times larger than the TD-DFT oscillator strengths. The next two states, D 3 and D 4 , which are again degenerate, are located at ∼1.8 eV at the TD-DFT level and at ∼ 2.0 eV at the semiempirical CISD level (these values hold for both geometries). The oscillator strengths for transitions to the D 3 ∕D 4 states are much larger than those for D 1 ∕D 2 . Again, the semiempirical oscillator strengths are larger than the TD-DFT ones, by a factor of ∼ 3. Starting from D 5 , there are differences in state ordering between TD-DFT and semiempirical CISD. Namely, D 5 and D 6 are degenerate at the semiempirical level, whereas at the TD-DFT level D 6 and D 7 are degenerate instead. However, the energy gap between D 5 and D 6 is very small, e.g. ∼ 0.01 eV for TD-DFT@DFT. The higher lying transitions (to D 8 , D 9 and D 10 ) are located at ∼3.8 eV at the TD-DFT level and have relatively large oscillator strengths. At the semiempirical CISD level, the D 0 → D 8 transition is the last transition with the excitation energy below 4 eV and this transition is rather intense. The D 9 and D 10 states appear at much larger excitation energies ≳ 6 eV because of the reduced active space. However, these states will not be used as initial states for the SH simulations (see below). The dominant orbital pairs contributing to the transitions are shown in Tables S2-S4 (for calculations at the DFT geometry). There, it is seen that (i) the character of D 1 -D 4 and D 8 states is very similar with both TD-DFT and the semiempirical CISD method, (ii) a state reordering for D 5  Table 2 Excitation energies and oscillator strengths for the lowest ten transitions of the adamantane cation calculated with ROOM3/CISD(9 × 10 ) and TD-UB3LYP/6-31G* at the UB3LYP/6-31G* and ROOM3/ CISD(9 × 10 ) optimized geometries Transition OM3/CISD(9 × 10 )@DFT TD-DFT@DFT OM3/CISD(9 × 10) @OM3/CISD(9 × 10) TD-DFT@OM3/ CISD(9 × 10)

Electronic spectra
Further, recently Crandall et al. reported the experimentally measured optical spectrum of Ada + obtained by EPD spectroscopy [39]. And even more recently, Kappe et al.
reported the electronic spectrum of cationic adamantane obtained by spectroscopy in helium droplets [40]. The characteristic feature of these spectra is the absence of vibronic progression at excitation energies smaller than 3.5 eV. At higher energies, on the contrary, Crandall et al. observed several well resolved peaks on top of the broad feature [39]. In order to compute vibrationally resolved absorption spectrum of Ada + , we first tried to do the excited state optimizations using TD-UB3LYP/6-31G*, but they failed. It is probably because linear response TD-DFT potential energy surfaces may be warped near a Jahn-Teller degeneracy [69]. In this respect, it should be pointed out that the spectra presented in this section do not include nonadiabatic effects (see also below). In passing we note that problems in excited state optimizations have also been encountered for polyaromatic hydrocarbons, and the definitive reason for this is not known [70]. Therefore, we resorted to a simple Vertical Gradient (VG) approximation [58] (which can be viewed as a version of what is known as IMDHO, Independent Mode Displaced Harmonic Oscillator -the model which does not include frequency alteration and Duschinsky rotation [71]). The vibrationally resolved absorption spectrum of Ada + in the VG approximation, calculated at the TD-UB3LYP/6-31G* level, is shown in Fig. 3. In qualitative agreement with experiments we observe two relatively intense bands, with maxima at ∼1.8 eV and ∼3.8 eV, respectively, and a lower energy weak band with a maximum at ∼ 1.1 eV. The lower energy band is made of the transitions to the D 1 and D 2 states (which are degenerate at the D 0 minimum geometry and remain degenerate in the VG approximation). The experimentally determined band maximum is at 1.12 eV [40], which agrees very well with the calculated band maximum position of 1.149 eV. The first intense band is made of the transitions to the D 3 and D 4 states (which are again degenerate). This band shows mostly unresolved vibronic structure, with the exception of a peak at the onset of the band, at 1.40 eV. This peak corresponds to the 0-0 transition between D 0 and D 3∕4 state (in the VG approximation). The second intense band is composed of transitions to D 8∕9 (which is degenerate) and to D 10 . Whereas the D 8∕9 band is structureless, the D 10 band shows some vibronic progression. It is in qualitative agreement with the spectrum reported by Crandall et al. [39]. Quantitatively, the first sharp peak of this band appears at 3.35 eV in the computed spectrum, whereas it is located at 3.60 eV in the experimental spectrum. This peak corresponds to the 0-0 transition between states D 0 and D 10 .
Furthermore, we computed the photoelectron (photoionization) spectrum in the VG approximation (Fig. 4). In addition, since we were able to optimize the cationic ground state ( D 0 ), using UB3LYP/6-31G*, we computed the S 0 → D 0 vibrationally resolved spectrum at the adiabatic Hessian (AH) level [58], also known as IMDHOFAD model [71], Independent Mode Displaced Harmonic Oscillator with Frequency Alteration and Duschinsky rotation (see the blue curve in Fig. 4).
In qualitative agreement with the earlier experiment [22] we observe three bands in the total PE spectrum (in the range below 14 eV), see grey line in Fig. 4. The first band, located at 9-10 eV, includes contributions from the D 0 , D 1 , and D 2 states. The second band (10-11.5 eV) comprises states D 3 -D 7 . And the third band (12-13.5 eV) includes D 8 -D 10 . The experimental band centers are 9.75, 11.25, and 13.40 eV, respectively [22]. We note that in contrast to the cationic absorption spectrum described above, the PE spectra of all For the S 0 → D 0 transition, we also computed the spectrum using the AH model. In comparison to the S 0 → D 0 VG spectrum, the AH spectrum is red-shifted and, moreover, the vibronic progression is different.
Finally, we note from the lowest band of the total spectrum, that the shift between 0-0 transitions of the S 0 → D 0 and S 0 → D 1∕2 is only ∼0.03 eV (compare the " D 0 " and " D 1 "/"D 2 " curves in Fig. 4). This small shift can already be anticipated from a single point TD-UDFT calculation for the cation at the S 0 minimum geometry (see Table S5), which gives the vertical excitation energy of the D 0 → D 1∕2 transition of ∼0.09 eV. A similar excitation energy of 0.07 eV computed with TD-B97X-D/6-311++G(2d,p) was reported by Candian et al. [45]. We should note, however, that the S 0 minimum geometry belongs to the T d symmetry, and the cationic ground state should be triply degenerate [72]. Therefore, the small spectral splitting should be the result of the unbalanced treatment of the ground and linearresponse states by linear response TD-DFT [73].
Moreover, we should stress that the spectra described here, besides the approximations discussed above, do not account for nonadiabatic effects, which are expected to be sizable for the adamantane cation, since it represents a Jahn-Teller problem [29].

Dynamics in D 0 with excess energy
We modeled adiabatic dynamics of the cation in the D 0 state with excess energy (see Fig. 1a,b). Specifically, we performed simulations with E h = 3, 4, and 10 eV for method (a) (excitation out of D 0 ) and E h = 17, 27, and 70 eV for method (b) (excitation out of S 0 ). For the latter, taking into account the vertical IP of 9-10 eV, the given excitation energies translate to 7-8, 17-18, and 60-61 eV, respectively, with respect to D 0 , providing the excess energy in the case of method (b). The excess energies of 3 and 4 eV for method (a) correspond to a higher energy part of the EPD spectra in [39]. For these energies, no fragmentation was observed in our simulations (see below) and the higher energy of 10 eV was tested in addition. For method (b), 70 eV corresponds to the nominal kinetic energy of electrons in typical EIMS experiments, and the energy of 27 eV is a rough estimate for the energy to be actually transferred to a molecule in these experiments (assuming internal excess energy per atom of 0.65 eV [61] and IP of 10 eV). We have also tried a lower energy of 17 eV to check that no fragmentation occurs when reducing energy in method (b).
The last snapshots from the MD trajectories, at time 10 ps (or earlier if a trajectory was aborted, which occurred to some trajectories at large excess energies) were analyzed for fragmentation. The relative abundance of all (neutral and cationic) fragments was calculated. For a fragment of a unique mass m we computed the relative abundance as the ratio of the number of fragments of mass m generated in the swarm of all trajectories to the total number of all fragments generated in this swarm. The relative abundance computed for various excess energies is shown in Fig. 5.
The left column of Fig. 5 corresponds to method (a) and the right column to method (b). For method (a), for E h = 3 and 4 eV, only the parent ion ( m = 136 u) is observed. For the larger energy of 10 eV, two trajectories out of 101 showed fragmentation. One of these two trajectories yielded H loss, and the other one showed fragmentation in C 4 H 8 and C 6 H 8 (see Fig. 5, bottom left). Moreover, 39 % of non-fragmented trajectories demonstrated rearrangement of the adamantane cage as exemplified by a derivative of the six-membered ring in Fig. 5, bottom left. We note that there are further types of the rearranged structures, e.g., the opened cage like the one shown in Fig. 5, top right. Interestingly, cage-opening of the amino derivative of Ada + , the amantadine cation (Ama + ), was recently investigated by George and Dopfer by means of infrared photodissociation Fig. 4 Photoelectron spectrum in the VG approximation (individual S 0 → D i contributions and their sum) as well as the S 0 → D 0 spectrum at the AH level. Calculations are done at the (TD-U)B3LYP/6-31G* level of theory spectroscopy and quantum chemical calculations [47,48]. It was found that amination largely lowers the barrier for the cage-opening reaction [48], from ∼1.9 eV for Ada + to ∼1.3 eV for Ama + (including H-shift, see [47]). However, given sufficient excess energies, we observe the cage-opening for Ada + in our simulations.
For method (b), fragmentation was not observed for E h = 17 eV (7-8 eV w.r.t. D 0 ), but 12 % of the parent ion signal at m = 136 u corresponds to the rearranged structures. For E h = 27 eV (17-18 eV w.r.t. D 0 ) and E h = 70 eV (60-61 eV w.r.t. D 0 ), there were quite a number of fragmentations. At 27 eV we got numerous fragments of C n H m , where n is ranging from 2 to 10, while at 70 eV the larger fragments were reduced to smaller fragments. Here the majority of the fragments were in the range of 1-3 carbon atoms per fragment. It agrees with a simple expectation that larger excess energy will induce more fragmentations and the fragments will be smaller in size.
We should also note that in methods (a) and (b) described above we completely neglected the kinetic energy of Ada or Ada + before excitation. Using arguments of thermodynamics, we can estimate this initial kinetic energy E ini kin from the initial temperature T ini (20 K for method (a) and 500 K for method (b)) as E ini kin = 3N−6 2 kT ini , where N is the number of atoms in the molecule ( N = 26 for Ada or Ada + ). For T ini = 20 K we then obtain E ini kin ≈ 0.06 eV, and for T ini = 500 K E ini kin ≈ 1.55 eV. The new excess energies can be obtained by subtracting E ini kin from the old excess energies, which, in turn, were defined as E h for method (a) and E h − IP for method (b). On the other hand, the initial conditions prepared by sampling the Wigner function (at 0 K) possess the kinetic energy E ini kin = 3.1 ± 0.7 eV. We performed calculations modifying the method (a) to use the scaled Wigner initial velocities instead of the randomly generated ones. Taking into account the initial kinetic energy of about 3.1 eV, the photon energy of 3 eV will result in an excess energy of about  Fig. 1a, b). Selected fragments are also shown 6.1 eV. The simulations performed with this excess energy showed no fragmentation within 10 ps also (see Fig. S1).
Overall, the MD simulations described above (which are performed in the D 0 state only) suggest that almost 10 eV of excess energy is required to fragment the adamantane cation within 10 ps.

Surface hopping simulations
In order to model internal conversion and explicitly account for the effect of cationic excited states, we performed nonadiabatic dynamics simulations for Ada + using Tully's surface hopping approach in combination with the ROOM3/CISD(9 × 10 ) electronic structure method. First, the electronic, vertical absorption spectra of the cation were calculated for the geometries sampled from the Wigner function. The stick spectra were broadened with unnormalized Gaussians as: Here, I is the intensity, E the excitation energy, E i, and f i, are the calculated excitation energy and oscillator strength, respectively, for the D 0 → D i transition of snapshot , is a broadening parameter ( = 0.05 eV in this work), and M is a normalization factor (here, M is chosen to be the maximum intensity of the broadened spectrum). The absorption spectrum of Ada + computed this way is shown in Fig. 6.
In the total, broadened spectrum, we can identify the following bands: (1) the lower energy weak band composed of the transitions to D 1 and D 2 , (2) the intense band including the transitions to D 3 -D 7 , (3) the band stemming from the D 0 → D 8 transition, and (4) the band corresponding to D 9 . We note that due to geometrical distortions of initial geometries, all transitions may acquire nonzero oscillator strengths. Comparing to Fig. 3, we see that both spectra exhibit two intense bands in the excitation energy range below 4.5 eV. We then modeled nonadiabatic dynamics after excitation to bands (1)-(3). As discussed above, the excitation energy of D 9 is too large when using the ROOM3/CISD(9 × 10 ) approach. Therefore, we do not consider dynamics initiated by excitation into band (4). Specifically, four batches of SH simulations were considered: (i) excitation to the D 1 ∕D 2 states, (ii) excitation to the D 3 ∕D 4 states, (iii) excitation to an energy window (2.3, 2.9) eV, and (iv) excitation to the D 8 state, as described already in the Methods section. The electronic state populations, calculated as fractions of trajectories in the respective state, are shown for all the cases in Fig. 7, for the first 200 fs (populations for the total propagation time of 1 ps are shown in Figs. S2-S5). We observe an ultrafast internal conversion to the ground state ( D 0 ). Using a single exponential function of the from P D 0 = 1 − exp(−t∕ ) to fit the ground-state population we obtain time constants of 10, 39, 37, and 100 fs for cases (i)-(iv), respectively. Thus, the low-lying excited states of Ada + exhibit ultrashort lifetimes. In particular, for case (iv) (excitation to D 8 ) using monoexponential fitting for D 8 , P D 8 = exp(−t∕ ) , we find the lifetime = 56 fs. Interestingly, Crandall et al. deduced the lifetime of 32 ± 11 fs from their spectrum [39], which is of the same order of magnitude. We also note that the ultrashort time constant of 10 fs when exciting to D 1 ∕D 2 is similar to the time scale of the D 1 → D 0 electronic relaxation in the ethylene cation ( ∼ 7 fs) [74].
Further, we analyzed the SH trajectories for fragmentation. Remarkably, all four batches of the SH simulations demonstrated the viability of the ultrafast H loss channel at the used (low) excitation energies, in contrast to the simulations in the D 0 state only. The H loss yield is, however, small, at least for the period of simulation (1 ps), see Fig. 8. There we plot the H loss yield, defined as the ratio of the number of trajectories showing hydrogen dissociation to the total number of trajectories, as a function of excitation energy using bars centered at the mean excitation energy and having the width of twice standard deviation (for the initial excitation energy). After excitation to D 1 and D 2 , the H loss yield is 0.9 %. For excitation to D 3 and D 4 , it is 1.5 %. Similarly, for excitation to the energy window (2.3, 2.9) eV, the yield is 1.8 %. After exciting D 8 , we find an H loss yield of 7.8 % (Fig. 8). All these numbers refer to the time interval of 1 ps of dynamics. Furthermore, we find that the H loss takes place in the ground state ( D 0 ), i.e. after internal conversion. At that, times of fragmentation span the whole interval of 1 ps. An example of the reactive trajectory (showing H dissociation) is presented in Fig. 9. It is seen that, for this particular trajectory, internal conversion to the ground state takes about 20 fs, whereas additional ∼ 20 fs are required for H loss to happen. To understand why H loss happens at low excitation energies in our SH simulations, we estimated the H dissociation energy in the D 0 state using ROOM3/CISD(9 × 10 ) and UOM3 (with Fermi smearing). The dissociation energy was calculated as the difference between the ground state minimum of C 10 H + 16 and a structure with an elongated by 10 Å CH distance ( r 1 in Fig. 2), which was optimized under the constraint of fixed r 1 . The ROOM3/CISD(9 × 10 ) dissociation energy was found to be 1.02 eV, whereas UOM3 (with Fermi smearing) yielded 2.09 eV. The CCSD(T) estimate is 1.12 eV [40]. In addition, we performed BOMD runs in the D 0 state, starting with an initial condition for one of the reactive SH trajectories (initiated in D 1 ) and using either ROOM3/CISD( 9 × 10 ) or UOM3+Fermi smearing. We observed H loss at the CISD  level but not at the UOM3 level. Thus, the low-energy H loss observed in our SH simulations can be explained by the lower ground-state H dissociation energy at the ROOM3/ CISD(9 × 10 ) level.
In addition, for the largest excitation energy (excitation to D 8 ), besides H loss, we observed also H 2 loss (only 2 trajectories) and C 3 H 5 loss (only 1 trajectory). We did not observe the formation of other larger fragments though (during 1 ps). However, we observe the rearrangement of the adamantane cage. Namely, we find 0.3 %, 0.5 %, 1.0 %, and 3.7 % of the parent signal coming from the rearranged structures for batches (i), (ii), (iii), and (iv), respectively. Finally, we note that in the breakdown diagram obtained by VUV iPEPICO spectroscopy [45], the parent ion signal is not present for energies > 11.65 eV with respect to S 0 (or > 2.36 eV with respect to D 0 ). Therefore, the corresponding fragmentation processes appear to be slower than the time scale addressed with the MD simulations in our work.

Summary and conclusions
We performed a computational study of the adamantane cation (Ada + ) addressing several aspects of its valence spectroscopy and molecular dynamics. First, we computed vibrationally resolved absorption and photoelectron spectra of Ada + and Ada, respectively, considering transitions not only to the lowest but also to the higher lying electronically excited states of the cation (Ada + ). These calculations were done at the (TD-)DFT (B3LYP/6-31G*) level of theory using the vertical gradient approximation. A semiquantitative agreement is found between the computed and experimental spectra. Second, we performed adiabatic ( D 0 ), ab initio molecular dynamics simulations employing the semiempirical OM3 electronic structure method and providing certain amount of the excess energy to simulate the fast fragmentation of Ada + after (photo)excitation. While the fast (on the time scale below 10 ps) fragmentation was observed for the excess energies of ≳ 10 eV with respect to D 0 , it was not detected for lower excess energies (at the UOM3+Fermi smearing level). Third, we applied the surface hopping approach to study internal conversion of Ada + . The nonadiabatic dynamics were modeled for 1 ps taking into account up to ten electronic states and using configuration interaction singles and doubles based on the OM3 molecular orbitals. We found that Ada + undergoes an ultrafast internal conversion to the D 0 state (on a time scale of 10-100 fs depending on initial excitation energy). The internal conversion was found to be followed by ultrafast fragmentation (predominantly H loss) for some trajectories. Specifically, we found the H loss yield of ∼1-8% for excitation energies ∼1-4 eV. The occurrence of the H loss in this case originates from the lower ground-state ( D 0 ) H dissociation energy at the ROOM3/CISD(9 × 10 ) level. Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.