Coronene: a model for ultrafast dynamics in graphene nanoflakes and PAHs

Abstract

Assuming a delta pulse excitation, quantum wavepackets are propagated on the excited state manifold in the energy range from 3.4–5.0 eV for coronene and 2.4–3.5 eV for circumcoronene to study the time evolution of the states as well as their lifetimes. The full-dimensional (102 and 210 degrees of freedom for coronene and circumcoronene respectively) non-adiabatic dynamics simulated with the ML-MCTDH method on twelve coupled singlet electronic states show that the different absorption spectra are only due to electronic delocalisation effects that change the excited state energies, but the structural dynamics in both compounds are identical. Breathing and tilting motions drive the decay dynamics of the electronic states away from the Frank–Condon region independently of the size of the aromatic system. This promising result allows the use of coronene as a model system for the dynamics of larger polycyclic aromatic hydrocarbons (PAHs) and graphene one dimensional sheets or nanoflakes.

---

1. Introduction

Photoinduced processes play a crucial role in many fields of chemistry, physics, biology, and material sciences, ranging from their involvement in vision to anticancer therapies and optoelectronics.^1–5 Understanding such processes is essential to be able to design new compounds and tune the optical and electronic properties of materials, enabling a desired functionalisation. Particularly in the field of optoelectronic devices such as solar cells, light-emitting diodes (LEDs), and lasers, photoinduced processes are fundamental, influencing their performance and efficiency. To conveniently tune their performance, one must understand the behavior of the excited states of the target materials and study their time evolution. Studying these processes enhances our understanding of fundamental physical phenomena, facilitates the development of new theoretical models, and guides the design of novel materials with tailored properties.^6–8

Polycyclic aromatic hydrocarbons (PAHs) are π-conjugated organic molecules that consist of multiple fused aromatic rings. They are naturally found as products from the incomplete combustion of organic materials, forest fires, and volcanic eruptions.^9,10 Leaving aside their threat to humans health,¹¹ in the last years they have become highly relevant in industrial applications such as the production of dyes and pigments and pharmaceuticals^12,13 due to their photostability and photophysical properties.^14,15 These compounds have also been found to be a key piece in the understanding the chemistry of the interstellar medium.¹⁶

Moving towards the solid state, graphene quantum dots (GQDs) are nanostructures composed of small fragments of graphene, which is a two-dimensional carbon allotrope, typically showing lateral dimensions ranging from a few to several hundred nanometers.¹⁷ GQDs possess unique optical, electronic, and chemical characteristics, distinct from bulk graphene, presumably arising from their size, edge structures, and quantum confinement effects.^18–20 In recent years, it has been proposed that their tunable photoluminescence may be related to the method of synthesis, attributed to doping effects, surface functional groups, and the presence of molecular fluorophores,^21–26 which has been widely investigated also by theoretical computations, recently reviewed by M. Langer et al. in ref. 27. As a consequence of their astonishing properties, GQDs have applications in optoelectronics,^22,28–30 energy storage,^31,32 sensing,^20,33 and biomedicine.^34,35

Graphene nanoflakes, also known as graphene nanoribbons or nanostrips, are narrow and elongated structures of graphene, and similarly to GQDs, they can be synthesized with different widths, lengths, and edge structures, resulting in diverse electronic properties.^36,37 Their edge-state-dependent properties are due to quantum confinement along the width of the nanoflakes. In general, they are larger than GQDs, and therefore the quantum confinement effects may be less pronounced.

Coronene and circumcoronene are both examples of polycyclic aromatic hydrocarbons that are often used as models for studying the properties of larger PAHs, GQDs, and graphene nanoflakes.³⁸ Coronene is a PAH consisting of a planar molecule with a hexagonal shape, composed of seven fused benzene rings. It has a well-defined structure and is often used as a model system to study the properties of extended aromatic compounds due to its relative simplicity. Coronene is known for its fluorescence properties and has been extensively studied for its electronic, optical, and chemical properties.³⁹ It is also used as a benchmark molecule in theoretical calculations and experimental studies.^40–43 UV-VIS results show a band at 4.1 eV, with a shoulder displaced to the blue (4.3 eV). A forbidden band, allowed only by symmetry breaking, appears at 3.6 eV.

Circumcoronene is a larger PAH composed of a circular arrangement of benzene rings. It is formed by the fusion of multiple coronene units, resulting in a larger and more complex structure. It contains exactly double the amount of carbon and hydrogen atoms than coronene (48 carbons and 24 hydrogens). Due to its size and extended π conjugation, circumcoronene exhibits unique electronic and optical properties. It is often employed as a model system to investigate the properties of graphene nanoflakes because of its structural similarity (and its larger size compared to coronene) and the possibility of tuning its electronic structure through chemical modifications. The synthesis of circumcoronene has always remained a challenge and it has been recently synthesized in solution,⁴⁴ although a UV-VIS spectrum is not yet available.

Both coronene and circumcoronene have been widely studied in various fields, including materials science, organic electronics, and nanotechnology, as they provide valuable insights into the behavior of PAHs and graphene nanoflakes. Due to their reduced size compared to the material structures of interest, one can use more accurate theoretical methods to deeply understand their properties. Therefore, they serve as important model systems for understanding the fundamental properties and potential applications of larger carbon-based materials.³⁸ Although there are excited state dynamical studies on smaller PAHs, such as benzene,^45,46 naphthalene dimer,⁴⁷ anthrazene⁴⁸ and pyrene,⁴⁹ the field is lacking an exhaustive work on the excited state behavior of coronene and circumcoronene. The purpose of this article is to compare the excited state dynamics on coronene and circumcoronene molecules, to assess the feasibility of using the coronene molecule as a model system for larger polycyclic aromatic hydrocarbons.

2. Theoretical and computational framework

2.1 Electronic structure details

To compute vertical energies, gradients and couplings needed throughout this work, we benchmarked several TD-DFT functionals and basis sets against the wave function method EOM-CCSD/cc-pVDZ for the coronene molecule. For the CCSD ground state reference, the symmetry was chosen to be D_2h, the highest symmetry available. The Q-Chem package⁵⁰ was used to compute the reference EOM-CCSD vertical energies, whereas and Gaussian 16⁵¹ was chosen to compute the vibrational frequencies of both systems. For every other calculation, i.e., Wigner-displaced geometries, parametrisation of the vibronic coupling Hamiltonian and on-the-fly dynamics, we used the ORCA 5.0 electronic structure program.^52,53 The details about the benchmark we performed are included in the ESI,† and they are supported by a previous work by B. Shi et al.³⁸ The TD-B3LYP-D3/def2-SVP method was proven to provide an excellent agreement with respect to experimental results while keeping the computational time feasible. For circumcoronene, we used the same combination to compute the parametrised potentials: TD-B3LYP-D3/def2-SVP.

2.2 The vibronic coupling model

To describe the potential energy surfaces of the systems, we constructed an approximate Hamiltonian as a first order Taylor expansion of the vertical electronic state energies calculated at the Frank–Condon (FC) point.

A diabatic basis ensures a smooth description of the potential energy surfaces and mass- and frequency-weighed normal mode coordinates Q_α(α = 1, …, f), where f is the number of vibrational degrees of freedom (DOFs) of the system, are chosen as the coordinate basis to expand the potential energy surfaces. Tuning modes carry the gradient, driving the molecules away from the FC region, whereas coupling modes mix pairs of electronic states and are responsible for the population transfer. When this expansion is truncated in the first order, the method is known as the linear vibronic coupling (LVC) model, and the Hamiltonian looks as follows:

H = H⁽⁰⁾ + W⁽⁰⁾ + W⁽¹⁾. (1)

The zeroth-order Hamiltonian H⁽⁰⁾ is defined by the diagonal kinetic energy operator of the N electronic states to be included,A simple way of defining the zeroth-order potential matrix W⁽⁰⁾ is to use the harmonic approximation shifted by the vertical excitation energies E_i at the Frank–Condon (FC) geometrywhile W⁽¹⁾, the first-order corrections of the potential matrix are defined using the following expansion:The κ⁽ⁱ⁾_α linear intra-state vibronic couplings carry the gradients (along the α tuning modes) of the i excited potential energy surface at the FC point. The λ^(ij)_α linear inter-state vibronic constants measure the coupling between each pair of electronic states i and j along the normal mode α.

In this work, we used the LVC model as implemented in the SHARC molecular dynamics package,^54–56 which generates all the LVC model parameters⁵⁷ from (a) the calculation of the ground-state frequencies at the gas-phase optimised geometries of coronene and circumcoronene and (b) 2(3N − 6) + 1 excited state energy calculations of displaced geometries in positive and negative directions along every vibrational mode from the equilibrium configuration. From these computations, diabatic state gradients and couplings were computed from numerical derivatives, as detailed in ref. 57. As we are dealing with systems of non-abelian symmetry groups, a rotation of the Hamiltonian was applied to the Hamiltonian to correct spurious couplings, as described in ref. 58.

2.3 The ML-MCTDH wavepacket propagation method

The nuclear dynamics of the system can be studied by solving the time-dependent Schrödinger equation (TDSE)where the Hamiltonian operator H can be generalized to describe a system characterized by a set of coupled potential energy surfaces in the diabatic picture, like in eqn (1). The total wavefunction of the system is a function of the normal modes and its definition depends on the method we use to solve this first-order differential equation.

A common approach to solve the TDSE is the multiconfigurational time dependent Hartree (MCTDH) method,^59,60 which expands the total wave function in a set of time-dependent basis functions ψ^(α)_p as

which are known as single-particle functions (SPFs) and are expanded in a time-independent primitive basis set χ^(α), constructed using a DVR⁶¹ grid,Due to the exponential scaling of the computational effort with respect to the dimensionality of the system, MCTDH is normally used for systems with no more than 30 vibrational degrees of freedom. For higher-dimensional systems, the multi-layer MCTDH (ML-MCTDH)⁶² approach can be used, where a recursive layering structure is constructed in a way that the SPFs are expanded in a set of time-dependent SPFs and only the last layer is expanded using the time-independent DVR primitive basis, still ensuring perfect convergence to the exact quantum dynamics result at the nuclear basis limit.

All these methodologies are the core of the Quantics package,^63,64 used throughout this work to study the quantum dynamical properties of coronene and circumcoronene.

2.4 Trajectory surface hopping

The Tully surface hopping (TSH)^65,66 method propagates the nuclei following the Newton equations of motionwhere the atom A is moving in a single potential energy surface (V_i(R(t))) pushed by forces corresponding to its gradient. Compared to MCTDH, instead of a full nuclear wavepacket, TSH consists of individual trajectories that move completely uncoupled from the bundle.

The propagation from one-time step to another is performed using the velocity Verlet algorithm,⁶⁷ where velocities v_A(t) and the atomic coordinates R_A(t) change according to

The nuclei are assumed to move adiabatically most of the time, and in some particular regions where the coupling with other electronic states is high, undergo non-adiabatic transitions known as hops, thus the name surface hopping.^65,66 At every time step, the nuclei have the choice to change the electronic state based on a probability to hop that depends on the coefficients of the electronic states involved. The electronic wavefunction is therefore expressed as a sum of different electronic states with time-dependent coefficients using the Born–Huang expansion. Inserting this wavefunction ansatz in the TDSE yields to the following equations for the propagation of the coefficients:where the Hamiltonian term is directly calculated from the precomputed PES or via electronic structure methods, and the second term is calculated from the NAC between electronic states and velocities: K_βα = 〈Ψ^el_β|∇_R|Ψ^el_α〉v_R. Throughout this work, we have used the SHARC implementation of TSH⁶⁸ both calculating electronic quantities on-the-fly and on precomputed potentials, as well as the more accurate wavepacket method ML-MCTDH on the same potentials. The reason for using the less-accurate method TSH was twofold: first, to assess the role of triplets since the spin orbit coupling was negligible at the Frank–Condon geometry and, second, to compare whether dynamics on precomputed potentials were a good approximation to the dynamics computed on-the-fly. The TSH calculations used the “fewest switches” algorithm,⁶⁶ along with the decoherence correction of Granucci and Persico.⁶⁹ The time-step was chosen to be 0.5 fs. The hopping probabilities were obtained from state overlaps.⁷⁰ 400 initial structures and velocities were generated from a Wigner distribution. More system-specific computational details are provided at the beginning of the corresponding sections.

3 Results: quantum dynamics of coronene

3.1 Vertical absorption and static results

Using the level of theory introduced in Section 2.1, we computed the corresponding vertical excitations from the singlet ground state (X, S₀) of coronene, collected in Table 1 in energetic order following the TD-B3LYP-D3 results. Besides the experimental ref. 41, Table 1 includes also the electronic character of every excited state, providing the molecular orbitals (MOs) with the largest contributions. The graphical representations of these MOs are included in Fig. 3, where the last 4 occupied orbitals are in the bottom figure, and the first 4 virtual orbitals are in the top row. The MOs represented in Fig. 3 are degenerate in the D_2h symmetry group in the following pairs of MOs: (LUMO+3, LUMO+2), (LUMO+1, LUMO), (HOMO, HOMO−1), (HOMO−2, HOMO−3).

Table 1 Vertical excitation energies, in eV, and their corresponding oscillator strengths obtained using EOM-CCSD/cc-pVDZ level of theory and TD-B3LYP-D3/def2-SVP for coronene in D_2h symmetry. The experimental absorption maxima are taken from ref. 41

State^a	Symm.^b	Character^c	TD-B3LYP-D3/def2-SVP	EOM-CCSD/cc-pVDZ	Expt.^41
a Adiabatic state labeling in energetic order at the FC region based on the TD-B3LYP-D3/def2-SVP results. b Diabatic symmetry labels within the D2h molecular symmetry group of the singlet excited states (at the TD-B3LYP-D3/def2-SVP level). c The transitions between MOs of coronene that contribute the most (>90%) to each excited state. The symmetry labels correspond to the D2h molecular symmetry group.
S0	X(Ag)		0.00	0.00
S1	1B2u		3.22 (0.00)	4.06 (0.00)
S2	1B3u		3.51 (0.00)	3.30 (0.00)	3.65 (0.51)
S3	1B1g		4.32 (0.00)	4.55 (0.00)
S4	1Ag		4.32 (0.00)	4.55 (0.00)
S5	2Ag		4.35 (0.00)	4.82 (0.00)
S6	2B2u		4.38 (0.85)	4.85 (1.14)	4.09 (2.24)
S7	2B3u		4.38 (0.85)	4.85 (1.14)	4.09 (2.24)
S8	2B1g		4.42 (0.00)	4.89 (0.00)
S9	3B1g		4.58 (0.00)	5.27 (0.00)
S10	4B1g		4.70 (0.00)	5.34 (0.00)
S11	3Ag		4.70 (0.00)	5.27 (0.00)

---

Among all the vertical transitions involving these orbitals, only two of them show oscillator strength at 4.4 eV, which corresponds to the excitation to the 2B_2u (S₆) and 2B_3u (S₇) singlet excited states, matching the maximum absorption peak of the experimental spectrum around 4 eV. These two bright states are degenerate in D_2h, and they correspond to a E_1u state in D_6h.³⁸ There is another absorption peak in the experimental spectrum at 3.6 eV which appears dark in the FC region and corresponds to the excitation to the 1B_3u (S₂) state. The transition to this state is forbidden by symmetry, but it becomes accessible through the vibronic couplings between states. This effect has been observed by generating a set of initial conditions displaced from the equilibrium configuration using a Wigner distribution based on the S₀ state frequencies computed at the TD-B3LYP-D3/def2-SVP level of theory. The vibrationally-resolved absorption spectrum obtained by exciting from all the initial configurations is represented in the upper panel (a) of Fig. 1, where a new bright band corresponding to excitation to the 1B_3u (S₂) state is observed around 3.4 eV. Only singlet states are included in this work due to the negligible Spin–Orbit Couplings (SOCs) present at the equilibrium configuration. The initial parametrisation of the Hamiltonian included also twenty triplet states, i.e., sixty electronic states more. Every matrix element was very small, indicating that the triplet states will not be populated at this geometry. To further analyse whether at different geometries changes that statement, fifty four TSH on-the-fly trajectories were propagated for 20 fs, in which no population transfer to any of the triplet excited states is observed, at least at the short time scales analysed in this study.

Fig. 1 Computed absorption spectrum of coronene. Top panel (a) shows the static spectra computed by vertical excitations from the FC region using TD-DFT (black line) and including vibronic coupling by computing vertical excitations from 200 initial conditions selected from a Wigner distribution. The experimental absorption spectrum is taken from ref. 41. Lower panel (b) represents the dynamic absorption spectra obtained from the Fourier Transformation (FT) 〈Ψ(t)|Ψ(t = 0)〉 of the autocorrelation function of the quantum dynamics starting from 1B_3u (S₂) excited state computed with ML-MCTDH and including the 102 normal modes of coronene.

After the static analysis of the excited states in coronene, an LVC Hamiltonian, as described in Section 2.2, has been constructed in order to approximate the analytical potential energy surfaces including twelve coupled singlet states projected onto the 102 nuclear degrees of freedom of the molecule.

3.2 Excited state quantum dynamics from the first absorption band: 1B_3u (S₂)

To study the dynamical evolution of the diabatic state manifold, excited state dynamics have been carried out with different methods from both bright excited states, S₂ and S₆/S₇, as well as from the dark states to show how well it compares to non-adiabatic dynamics in the circumcoronene system.

For the case of the TSH, the 54 on-the-fly trajectories were initialised based on oscillator strengths using an energetic window of 3.0 to 3.75 eV to excite to the first “forbidden” band and performed on four coupled singlet electronic states (S₀–S₄) and six triplets (T₁–T₆) to minimise the computational effort needed. The TSH dynamics on LVC potentials started from the same geometries as the on-the-fly dynamics and included 54 trajectories on twelve singlet and eleven triplet electronic states. The ML-MCTDH dynamics were initialised as the Hartree product of the vibrational eigenstate of every degree of freedom. The dynamics on the twelve coupled electronic states was propagated using Runge–Kutta 5 for the nuclear degrees of freedom and short iterative Lanczos (SIL) for the electronic one. 4.6 × 10⁷ primitive basis functions were used for the grid and the calculations included 1.2 × 10³⁵ configurations. For the 8D reduced MCTDH dynamics the Bulirsch–Stoer extrapolation scheme was used for the propagation of the SPFs and the SIL integrator for the A-coefficients and the electronic propagation. The total size of the grid for the converged results in this system was 3.25 × 10⁹ and used 789 nuclear configurations.

The experimental gas phase absorption spectrum of coronene was recorded by S. Hirayama et al.,⁴¹ where the first absorption band of coronene was reported around ∼3.6 eV, which corresponds to excite the system to the second singlet excited state, 1B_3u (S₂). The population transfer from this excited state is shown in Fig. 2 using the three different methods for propagating the nuclear dynamics: (a) on-the-fly TSH classical trajectories, (b) TSH classical trajectories on the LVC potentials and (c) ML-MCTDH quantum dynamics on the same LVC potentials.

Fig. 2 Diabatic population transfer from the 1B_3u (S₂) excited state of coronene using three different propagations: (a) On-the-fly Tully Surface Hopping classical trajectories (TSH). The diabatisation has been made just tracing back to the state character of the initial populated state and populations are collected as the square of the quantum amplitudes of the states, (b) TSH classical trajectories using the LVC model for the PES and, (c) ML-MCTDH using also the LVC model potentials. Insets of molecular structures of coronene visited during both the on-the-fly and on LVC potentials using the TSH method are included in the figures (a) and (b) with 8892 and 2213 overlapped structures respectively.

In Fig. 2, the main difference between panels (a) and (b) is the PESs on which the classical nuclear trajectories were propagated. The more similar those panels, the better the LVC approximation in the PES construction. It is important to note that diabatic populations in panel (a) are not truly diabatic, but that the states have been labelled according to the diabatic state distribution that was present at the initial geometries, in a diabatic by ansatz manner. To clarify this fact, the combination of the S₁ and S₂ adiabatic populations is included in both (a) and (b) panels of Fig. 2, were a good agreement can be found. TSH in this case does not give the same diabatic results on-the-fly or on precomputed potentials and both panels are substantially different. This may be surprising since the LVC model has been proven to be a good approximation in rigid systems with impeded torsions, as it is the case of the coronene molecule. The performance of the LVC approximation has been also tested by inspecting the molecular configurations visited by the nuclear dynamics, represented as overimposed xyz structures in Fig. 2 right hand side. This figure shows a fairly harmonic behavior of the coronene geometry during the classical trajectories of the nuclei both over the LVC potential and the ‘real’ ab initio potential. Another proof-of-principle of the LVC approximation in this system is the fact that the absorption spectrum obtained with the PES is in good agreement with respect to the experimental, represented in the lower panel (b) of Fig. 1. As we can see in Table 1, the MOs involved in S₁ and S₂ are the same (with different relative weights), and they are degenerated in D_2h. To unravel why the diabatic population transfer differs, we checked the orbital character at the geometries visited along the on-the-fly dynamics and concluded that the differences arise due to a trivial rotation of the degenerated MOs.

In Fig. 2, the difference between panels (b) and (c) lies in the quantum treatment of the atomic nuclei in MCTDH in contrast to classical propagations during the TSH/LVC approach. From the comparison of these two panels, we can conclude that the TSH/LVC approach is a good approximation that covers qualitatively the population changes. However, the fine print of the dynamics of coronene is only provided by the most accurate method ML-MCTDH.

For the ML-MCTDH method, all 102 vibrational modes of the system have been included. The ML-tree structure, available on the ESI,† has been optimized by looking at the largest frequency-weighted κ⁽ⁱ⁾_α gradients and λ^(ij)_α couplings among states which involved the initial excited state 1B_3u (S₂).⁷¹

By using the same arguments, we can analyze those vibrations (normal modes) that contribute the most to the structure relaxation and drive the nuclear dynamics.

The identification of the leading vibrations in the dynamics is crucial, not only for unraveling the nuclear dynamics of coronene but for validating this methodology for future work with larger systems like graphene nanoflakes, and to understand the absorption and emission properties of these systems in general. For the case of coronene, the decay from S₂ is mostly due to the totally symmetric breathing mode ν₂₄, which shows a negligible diabatic coupling between S₂ and S₁ but large gradient on the S₂ state. Another leading mode in this dynamics is ν₁₉, the scissoring mode exhibiting the largest diabatic coupling. Two cuts of the diabatic PES along these two modes, along with their corresponding diabatic couplings are represented in the top panels of Fig. 5. The full-dimensional (full-D) dynamics have been compared with several reduced-dimensionality models, where only the most important modes have been included in a MCTDH calculation. We found the most relevant 8 modes which are the main responsible for the spectral features of the main band and during the first 10 fs of the dynamics. We systematically added modes to this 8D model trying to find a larger system to recover the full-D dynamics up to 50 fs. This was not possible, since there are many modes which contribute with small couplings to the Hamiltonian and cannot be left out. The spectrum resulting from the 8D dynamics and state populations compared to the full D results are depicted in the SI. The eight modes included are the three main breathing modes (full breathing v24, inner benzene breathings v60 and v77), four tilting modes (v19, v64, v79 and v88) and a CH stretching mode (v108). Note that the numeration of the modes started at v7 after removing the six translations and rotations.

After the vertical excitation from the ground state (GS) to the 1B_3u (S₂) state, the diabatic population transfer from this excited state is represented in Fig. 6 (top left panel) to up to 50 fs, where it can be seen that the population is equally distributed between S₂ and S₁ states already after the first 15 fs. When classical trajectories are used for the nuclei motion, a population transfer is also observed from the S₂ to S₁ state. However, the time scales observed in panels (b) and (c) of Fig. 2 are slightly different. Slower time scales are often observed when running classical trajectories, since they need to move towards the intersection, whereas a wavepacket has grid functions i.e., the possibility of population, in every part of the conformational space. From this comparison, it is suggested that the quantum effects play an important role in the nuclear motion over the PESs. We could also extract that electronic populations may not be the best description when dealing with high symmetry molecules, since arbitrary rotations impede a correct assignment of state characters.

3.3 Excited state quantum dynamics from every electronic state

To complete the dynamical study while considering the results obtained from the analysis of the quantum dynamics from the 1B_3u (S₂), ML-MCTDH has been carried out to unravel the quantum dynamics from every electronic excited state of coronene in the relevant energy range (below 5 eV). By exciting the system to each of the first 8 singlet states, we followed the corresponding diabatic populations. The corresponding figures of these population transfers are collected in the ESI† (Fig. S9–S16).

Since the excited states 2B_3u (S₆) and 2B_2u (S₇) are responsible of the most intensive peak in the experimental absorption spectrum, the quantum dynamics from these states are especially relevant. Due to the symmetry reduction in this calculations, 2B_3u and 2B_2u, which are degenerate in D_2h, correspond to the same state in the real symmetry group D_6h. This is confirmed by the fact that the quantum dynamics from both excited states are almost identical. After the excitation to these states, ∼70% of the population is transferred during the first 50 fs. From 2B_3u (S₆), the decay seems slightly faster as compared to starting the dynamics on the 2B_2u (S₇). Similar behavior is found for the diabatic population transfer from 1B_1g (S₃) and 1A_g (S₄) states, which are also degenerated in D_2h. The likeness between them is stronger during the first 10 fs, going from S₃ and S₄ to the degenerated pair S₁₀ and S₁₁ respectively, caused by the energy proximity and coupling among these states after vibrational relaxation. The first two excited states are also populated but the decay is faster from S₃. The diabatic transfer from 2 A_g (S₅) and 2 B_1g (S₈) are the fastest, where more than 70% of the population is lost before 10 fs.

4 Results: excited state quantum dynamics of circumcoronene and comparison to coronene

If we think about coronene as the smallest portion of a graphene sheet including all the symmetry operations (D_6h), then circumcoronene is the next step towards the generalization of graphene retaining the highest symmetry. Thus, we extended our study to this system, with twice as many atoms (72) as coronene and 210 vibrational modes.

4.1 Vertical absorption and static results

For consistency, the same methodology and level of theory introduced in Section 2.1, and used for coronene in Section 3.1, has been used to characterize the diabatic PES of circumcoronene. The vertical excitations have been computed also for this bigger graphene nanoflake using TD-DFT and included in Table 2, similarly to the analysis of coronene in Table 2, except for the absence of the EOM-CCSD and experimental references. Although this system has been synthesized recently by R. Gershoni-Poranne et al.,⁷² there is not an experimental absorption and emission spectrum, yet. One very interesting outcome of this work is the fact that the molecular orbitals involved in the excitations are equivalent to those in coronene, however, the excited state energies lay one electronvolt beneath due to delocalization effects. Table 1 includes also the relation of each excited state with its counterpart in the coronene molecule. This analysis was done by looking at the MOs involved in every transition and comparing those orbitals with the ones involved in the coronene excitations. The graphical representation of these MOs is included in Fig. 3, where the four highest occupied and four lowest virtual MOs of coronene and circumcoronene have been plotted. These MOs in circumcoronene are pairwise degenerated, similar to the coronene orbitals. The symmetry of each molecular orbital has been deluded by comparison with respect to the equivalent MO in coronene. The correlation between the MOs in the two systems has been represented by crossing lines in Fig. 3. It is interesting that the MOs that contribute the most to the bright excited states in coronene, 2B_2u (S₆) and 2B_3u (S₇), correspond to the same MOs that are involved in the two excited states in cirumcoronene showing a non-zero oscillator strength, 2B_2u (S₃) and 2B_3u (S₄).

Table 2 Vertical excitation energies, in eV, and their corresponding oscillator strengths were obtained using TD-B3LYP-D3/def2-SVP level of theory for circumcoronene in D_2h symmetry

State	TD-B3LYP-D3/def2-SVP	Coronene (Table 1)
S0	0.00	X(Ag) (S0)
S1	2.24 (0.00)	1B2u (S1)
S2	2.43 (0.00)	1B3u (S2)
S3	3.09 (1.27)	2B3u (S7)
S4	3.09 (1.27)	2B2u (S6)
S5	3.10 (0.00)	1Ag (S4)
S6	3.10 (0.00)	1B1g (S3)
S7	3.13 (0.00)	2B1g (S8)
S8	3.19 (0.00)	2Ag (S5)
S9	3.21 (0.00)	3B1g (S9)
S10	3.31 (0.00)	3Ag (S11)
S11	3.31 (0.00)	4B1g (S10)

---

Fig. 3 Top view of molecular orbitals of circumcoronone showed in comparison with respect to the corresponding ones in coronene. All these orbitals correspond to π* orbitals asymmetric with respect to the molecular plane. Symmetry labels are included for the coronene system, connected with lines with equivalent MO in circumcoronene, where the symmetry can be extrapolated.

The circumcoronene absorption spectrum, represented in Fig. 4, has been simulated in the same manner we did for coronene in Fig. 1, where the upper panel shows the static spectrum with the vertical absorption peak allowed by symmetry (S₃ and S₄) and including then vibronic effects by computing vertical transitions from 200 initial configurations generated from a Wigner distribution, which reveals a smaller absorption peak at lower energies, corresponding to excite the system to S₂.

Fig. 4 Computed absorption spectrum of circumcoronene. Top panel (a) shows the static spectra computed by vertical excitations from the ground state optimised geometry using TD-DFT (black line) and including vibronic effects by computing vertical excitations from 200 initial conditions selected from a Wigner distribution, which are shown in colors after convoluting with gaussians with a width of 0.2 eV. Lower panel (b) represents the dynamic absorption spectra obtained from the Fourier Transformation (FT) of the autocorrelation function 〈Ψ(t)|Ψ(t = 0)〉 of the quantum dynamics starting from 1B_3u (S₂) excited state computed with ML-MCTDH and including the 210 normal modes.

4.2 Excited state quantum dynamics from every electronic state in the coupled manifold

The same procedure has been employed to study the out-of-equilibrium dynamics of circumcoronene to observe the effect of the enlargement of the system in the molecular plane.

By exciting the system to the second singlet excited state (1B_3u, S₂) the dynamic absorption spectrum has been recovered from the Fourier Transformation of the autocorrelation function and plotted in the bottom panel of Fig. 4 to compare better with respect to the static absorption. The quantum dynamics for this system have been studied using ML-MCTDH on the LVC precomputed potentials calculated at the TD-B3LYP-D3/def2-SVP including all the vibrational degrees of freedom. During these dynamics, the grid sizes were 5.91 × 10²¹⁹ and 2.43 × 10⁵⁰ nuclear configurations were included. Due to its large size, on-the-fly TSH dynamics is out of reach for this system.

We also performed a vibrational mode analysis by exploring the analogous modes that were relevant in coronene. The scissoring mode in coronene ν₁₉ is equivalent to ν₃₄ in circumcoronene, which is equally important for the de-excitation dynamics from S₂ in this molecule. In the same way, the totally symmetric breathing ν₂₄ in coronene is equivalent to ν₂₈ in circumcoronene, which is also responsible for leading the molecules away from the Frank–Condon region on the S₂ state. In Fig. 5, these four cuts of the adiabatic PES are represented in order to show the clear correlation between these important vibrations in both systems. Although the behavior of both PES is quite similar, it is important to emphasize here that they are displaced to lower energies and that there is a higher density of states in the case of circumcoronene, which can be explained by the size difference and the larger electronic delocalization in this system.

Fig. 5 Potential energy cuts along the two leading modes in the dynamics from the excited state 1B_3u (S₂) of coronene (top panels) and the corresponding modes in circumcoronene (bottom panels). The right panels correspond to the most important symmetric breathing modes with negligible coupling and large gradients, while the left panels correspond to the most relevant scissoring mode in both systems, exhibiting large diabatic couplings. The diabatic coupling between 1B_2u (S₁) and 1B_3u (S₂) is plotted in grey.

In Fig. 6 we grouped the excited state dynamics in coronene and circumcoronene from S₂ to S₈ according to their similarities. 1B_2u and 1B_3u states behave in a very similar manner in coronene and circumcoronene, transforming population among them in the sub 20 fs timescale. After this time, an equilibrium is reached and both states keep half the population. 1B_1g, 1A_g, and 2B_1g states have an equivalent initial decay lead by asymmetric in-plane motions. The population of the states is kept constant after 20 fs. When exciting to 2A_g in coronene and circumcoronene (S₅ and S₈ respectively) a complete population decay to other states is reached within 10 fs. The bright 2B_2u and 2B_3u states, however, decay at a slower pace, showing lifetimes above 250 fs.

Fig. 6 Population transfer from ML-MCTDH full-D dynamics starting in every excited state for coronene (left column) and circumcoronene (right column). The order of states has been altered to keep the same order as in Table 1 so that the dynamics are comparable.

5 Conclusions

We have benchmarked the level of theory to describe the electronic potential energy surfaces in coronene, finding that TD-B3LYP-D3/def2-SVP provides a good agreement with experimental results. A parametrised LVC model explains the symmetry-forbidden vibronic band at around 3.4 eV. The full dimensional excited state dynamics is driven by breathing and tilting modes. No reduced mode combination, however, leads to the full dimensional result, indicating that most of modes contribute to the dynamics with small couplings. On-the-fly excited state dynamics where electronic quantities are computed along the trajectories show that the contribution of triplet states is negligible and indicates that the conformational molecular changes are correctly described by displaced harmonic potentials. As the D_2h assumed symmetry mixes electronic states that would be degenerated at the FC region in the D_6h group, differences in population transfer arise in on-the-fly dynamics starting in the S₂ state when compared to precomputed potentials. The same LVC/TD-B3LYP-D3/def2-SVP methodology was applied for the circumcoronene molecule with 210 nuclear degrees of freedom. Although the absorption spectrum lays 1 eV below the spectrum of coronene, the spectral features are very similar. In fact, the excited state dynamics started from S₂ to S₈ involving twelve coupled electronic states and 210 DOFs agree very well with the coronene dynamics. This result shows that only one electronic structure calculation at the FC in the larger system and a state mapping to the smallest unit – coronene – are needed to describe the early ultrafast dynamics of polyaromatic hydrocarbons and their derivatives.

Data availability statement

The benchmark for the vertical energies in coronene, the comparison from spectra computed from 8D dynamics and full dimensionality for coronene, the ML-MCTDH trees which show the nuclear base function distribution among the 102 and 210 degrees of freedom and the detail from the dynamics from every electronic state for both molecules can be found in the ESI.† The frequency calculations at the optimised geometries of both compounds, quantics inputs, operator files and initial conditions for the on-the-fly TSH dynamics are in the comp data.zip file. The computational data has also been uploaded to the GREDOS repository (Universidad de Salamanca) and can be found at https://hdl.handle.net/10366/153091.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

SG acknowledges NextGenerationEU funds (María Zambrano Grant for the attraction of international talent) and the USAL grant “Programa Propio C1”. The authors thank the funding by Spanish Ministry of Science and Innovation (MCIN/AEI/10.13039/501100011033) grant no. PID2020-113147GA-I00.

Notes and references

1. W. Domcke and D. R. Yarkony, Ann. Rev. Phys. Chem., 2012, 63, 325–352.
2. M. Lan, S. Zhao, W. Liu, C.-S. Lee, W. Zhang and P. Wang, Adv. Healthcare Mater., 2019, 8, 1900132.
3. O. Ostroverkhova, Chem. Rev., 2016, 116, 13279–13412.
4. J.-Y. Xu, X. Tong, P. Yu, G. E. Wenya, T. McGrath, M. J. Fong, J. Wu and Z. M. Wang, Adv. Sci., 2018, 5, 1800221.
5. N. Keller and T. Bein, Chem. Soc. Rev., 2021, 50, 1813–1845.
6. T. C. Pham, V.-N. Nguyen, Y. Choi, S. Lee and J. Yoon, Chem. Rev., 2021, 121, 13454–13619.
7. E. Bassan, A. Gualandi, P. G. Cozzi and P. Ceroni, Chem. Sci., 2021, 12, 6607–6628.
8. J. M. Marin-Beloqui, S. Gómez, H. I. Gonev, M. Comí, M. Al-Hashimi and T. M. Clarke, Chem. Sci., 2023, 14, 812–821.
9. W. Giger and M. Blumer, Anal. Chem., 1974, 46, 1663–1671.
10. C. E. Cerniglia, Curr. Opin. Biotechnol, 1993, 4, 331–338.
11. K.-H. Kim, S. A. Jahan, E. Kabir and R. J. Brown, Environ. Int., 2013, 60, 71–80.
12. E. Huberman, H. Yamasaki and L. Sachs, Int. J. Cancer, 1976, 18, 76–82.
13. J. Berenbeim, S. Boldissar, S. Owens, M. Haggmark, G. Gate, F. Siouri, T. Cohen, M. F. Rode, C. S. Patterson and M. de Vries, Sci. Adv., 2019, 5, eaaw5227.
14. X. Shao, A. J. Aquino, M. Otyepka, D. Nachtigallová and H. Lischka, Phys. Chem. Chem. Phys., 2020, 22, 22003–22015.
15. F. Huisken, G. Rouillé, M. Steglich, Y. Carpentier, C. Jäger and T. Henning, Proc. Int. Astron. Union, 2013, 9, 265–275.
16. G. Malloci, G. Mulas and C. Joblin, Astron. Astrophys., 2004, 426, 105–117.
17. M. Bacon, S. J. Bradley and T. Nann, Part. Part. Syst. Charact., 2014, 31, 415–428.
18. J. Shen, Y. Zhu, X. Yang and C. Li, Chem. Commun., 2012, 48, 3686–3699.
19. S. Zhu, S. Tang, J. Zhang and B. Yang, Chem. Commun., 2012, 48, 4527–4539.
20. H. Sun, L. Wu, W. Wei and X. Qu, Mater. Today, 2013, 16, 433–442.
21. S. H. Jin, D. H. Kim, G. H. Jun, S. H. Hong and S. Jeon, ACS Nano, 2013, 7, 1239–1245.
22. T. Yuan, T. Meng, P. He, Y. Shi, Y. Li, X. Li, L. Fan and S. Yang, J. Mater. Chem. C, 2019, 7, 6820–6835.
23. K. Hola, Y. Zhang, Y. Wang, E. P. Giannelis, R. Zboril and A. L. Rogach, Nano Today, 2014, 9, 590–603.
24. J. Li, B. Wang, H. Zhang and J. Yu, Small, 2019, 15, 1805504.
25. O. Kozak, M. Sudolska, G. Pramanik, P. Cigler, M. Otyepka and R. Zboril, Chem. Mater., 2016, 28, 4085–4128.
26. M. L. Liu, B. B. Chen, C. M. Li and C. Z. Huang, Green Chem., 2019, 21, 449–471.
27. M. Langer, M. Paloncýová, M. Medved, M. Pykal, D. Nachtigallová, B. Shi, A. J. Aquino, H. Lischka and M. Otyepka, Appl. Mater. Today, 2021, 22, 100924.
28. F. Arcudi, V. Strauss, L. orević, A. Cadranel, D. M. Guldi and M. Prato, Angew. Chem., Int. Ed., 2017, 129, 12265–12269.
29. H. Tetsuka, A. Nagoya, T. Fukusumi and T. Matsui, Adv. Mater., 2016, 28, 4632–4638.
30. M. Semeniuk, Z. Yi, V. Poursorkhabi, J. Tjong, S. Jaffer, Z.-H. Lu and M. Sain, ACS Nano, 2019, 13, 6224–6255.
31. Q. Liu, J. Sun, K. Gao, N. Chen, X. Sun, D. Ti, C. Bai, R. Cui and L. Qu, Mater. Chem. Front., 2020, 4, 421–436.
32. X. Li, M. Rui, J. Song, Z. Shen and H. Zeng, Adv. Funct. Mater., 2015, 25, 4929–4947.
33. Z. Fan, S. Li, F. Yuan and L. Fan, RSC Adv., 2015, 5, 19773–19789.
34. T. H. Kim, J. P. Sirdaarta, Q. Zhang, E. Eftekhari, J. John, D. Kennedy, I. E. Cock and Q. Li, Nano Res., 2018, 11, 2204–2216.
35. S. Kadian, S. K. Sethi and G. Manik, Mater. Chem. Front., 2021, 5, 627–658.
36. A. Kuc, T. Heine and G. Seifert, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 085430.
37. C. Mansilla Wettstein, F. P. Bonafé, M. B. Oviedo and C. G. Sánchez, J. Chem. Phys., 2016, 144, 224305.
38. B. Shi, D. Nachtigallová, A. J. Aquino, F. B. Machado and H. Lischka, J. Chem. Phys., 2019, 150, 124302.
39. G. Bermudez and I. Y. Chan, J. Phys. Chem., 1986, 90, 5029–5034.
40. J. W. Patterson, J. Am. Chem. Soc., 1942, 64, 1485–1486.
41. S. Hirayama, H. Sakai, Y. Araki, M. Tanaka, M. Imakawa, T. Wada, T. Takenobu and T. Hasobe, Chem. – Eur. J., 2014, 20, 9081–9093.
42. T. E. Timofeeva, M. N. Egorova and A. E. Tomskaya, AIP Conf. Proceed, 2021, 2328, 050022.
43. D. Moran, F. Stahl, H. F. Bettinger, H. F. Schaefer and P. V. R. Schleyer, J. Am. Chem. Soc., 2003, 125, 6746–6752.
44. Y. Zou, X. Hou, H. Wei, J. Shao, Q. Jiang, L. Ren and J. Wu, Angew. Chem., Int. Ed., 2023, 62, e202301041.
45. D. Parker, R. Minns, T. Penfold, G. Worth and H. Fielding, Chem. Phys. Lett., 2009, 469, 43–47.
46. M. Galbraith, S. Scheit, N. Golubev, G. Reitsma, N. Zhavoronkov, V. Despré, F. Lepine, A. Kuleff, M. Vrakking and O. Kornilov, et al. , Nat. Commun., 2017, 8, 1018.
47. M. Gudem and M. Kowalewski, Chem. – Eur. J., 2022, 28, e202200781.
48. A. Marciniak, V. Despré, T. Barillot, A. Rouzée, M. Galbraith, J. Klei, C.-H. Yang, C. Smeenk, V. Loriot and S. N. Reddy, et al. , Nat. Commun., 2015, 6, 7909.
49. G. Braun, I. Borges, A. J. Aquino, H. Lischka, F. Plasser, S. A. Do Monte, E. Ventura, S. Mukherjee and M. Barbatti, J. Chem. Phys., 2022, 157, 154305.
50. Y. Shao, Z. Gan, E. Epifanovsky, A. T. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kus, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. DiStasio, H. Do, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. Hanson-Heine, P. H. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. D. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.-P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, E. Neuscamman, C. M. Oana, R. Olivares-Amaya, D. P. ONeill, J. A. Parkhill, T. M. Perrine, R. Peverati, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, S. M. Sharada, S. Sharma, D. W. Small, A. Sodt, T. Stein, D. Stück, Y.-C. Su, A. J. Thom, T. Tsuchimochi, V. Vanovschi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, J. Yang, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhao, B. R. Brooks, G. K. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. V. Voorhis, J. M. Herbert, A. I. Krylov, P. M. Gill and M. Head-Gordon, Mol. Phys., 2015, 113, 184–215.
51. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision C.01, Gaussian Inc., Wallingford CT, 2016.
52. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
53. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12, e1606.
54. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola and L. González, J. Chem. Theory Comput., 2011, 7, 1253–1258.
55. S. Mai, P. Marquetand and L. González, Int. J. Quant. Chem., 2015, 115, 1215–1231.
56. S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370.
57. F. Plasser, S. Gómez, M. F. Menger, S. Mai and L. González, Phys. Chem. Chem. Phys., 2019, 21, 57–69.
58. S. Gómez, M. Heindl, A. Szabadi and L. González, J. Phys. Chem. A, 2019, 123, 8321–8332.
59. M. H. Beck, A. Jäckle, G. A. Worth and H.-D. Meyer, Phys. Rep., 2000, 324, 1–105.
60. H.-D. Meyer, F. Gatti and G. A. Worth, Multidimensional Quantum Dynamics: MCTDH Theory and Applications, Wiley-VCH, Weinheim, Germany, 2009.
61. J. C. Light, I. P. Hamilton and J. V. Lill, J. Chem. Phys., 1985, 82, 1400–1409.
62. H. Wang and M. Thoss, J. Chem. Phys., 2003, 119, 1289–1299.
63. G. A. Worth, K. Giri, G. W. Richings, I. Burghardt, M. H. Beck, A. Jäckle and H.-D. Meyer, The QUANTICS Package, Version 1.1, University of Birmingham, Birmingham, UK, 2015.
64. G. Worth, Comput. Phys. Commun., 2020, 248, 107040.
65. J. C. Tully and R. K. Preston, J. Chem. Phys., 1971, 55, 562–572.
66. J. C. Tully, J. Chem. Phys., 1990, 93, 1061.
67. L. Verlet, Phys. Rev., 1967, 159, 98.
68. S. Mai, P. Marquetand and L. González, Int. J. Quant. Chem., 2015, 115, 1215–1231.
69. G. Granucci, M. Persico and A. Zoccante, J. Chem. Phys., 2010, 133, 134111.
70. F. Plasser, M. Ruckenbauer, S. Mai, M. Oppel, P. Marquetand and L. González, J. Chem. Theory Comput., 2016, 12, 1207–1219.
71. A. Lehr, S. Gómez, M. A. Parkes and G. A. Worth, Phys. Chem. Chem. Phys., 2020, 22, 25272–25283.
72. R. Gershoni-Poranne and A. Tsybizova, Angew. Chem., Int. Ed., 2023, e202305289.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp03656a

---