# A Density Functional Theory and Semiempirical Framework for Trajectory Surface Hopping on Extended Systems

Jan-Robert Vogt, Michael Schulz, Rafael Souza Mattos, Mario Barbatti, Maurizio Persico, Giovanni Granucci, Jürg Hutter, Anna Hehn

► To cite this version:

Jan-Robert Vogt, Michael Schulz, Rafael Souza Mattos, Mario Barbatti, Maurizio Persico, et al. A Density Functional Theory and Semiempirical Framework for Trajectory Surface Hopping on Extended Systems. Journal of Chemical Theory and Computation, 2025, 10.1021/acs.jctc.5c01082. [hal-05321933]

HAL Id: hal-05321933

https://hal.science/hal-05321933v1

Submitted on 20 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A Density Functional Theory and Semiempirical Framework for Trajectory Surface Hopping on Extended Systems

Jan-Robert Vogt, Michael Schulz, Rafael Souza Mattos, Mario Barbatti, Maurizio Persico, Giovanni Granucci, Jürg Hutter, and Anna Hehn*

<figure><img src="image_3.png" /></figure>

Cite This: https://doi.org/10.1021/acs.jctc.5c01082

<figure><img src="image_4.png" /></figure>

Read Online

## ACCESS

<figure><img src="image_5.png" /></figure>

Metrics & More

<figure><img src="image_6.png" /></figure>

Article Recommendations

<figure><img src="image_7.png" /></figure>

Supporting Information

**ABSTRACT:** Nonadiabatic molecular dynamics simulations provide a theoretical understanding of various excited-state processes in photochemistry, offering access to band widths, radiative or nonradiative relaxation and corresponding lifetimes, excited-state energies, and charge transfer. The range of method developments within the framework of time-dependent density functional theory is exceedingly large for molecular quantum chemistry. Still, it shrinks significantly when aiming to treat periodic boundary conditions. To address this gap and complement existing software packages for solid-state nonadiabatic molecular dynamics, we present an interface between the CP2K electronic structure and the NEWTON-X surface hopping codes. The interface features the generation of initial conditions, as well as adiabatic and non-adiabatic molecular dynamics, based on phenomenological or numerical time-derivative couplings. Setups are validated on gas-phase pyrazine, with electronic absorption spectra and excited-state populations for transitions between the lowest singlet states being in agreement with established molecular quantum chemistry methods. Extending the system size to crystalline pyrazine, limitations of approximate couplings are discussed, and the efficiency and applicability of the interface are demonstrated by computing broad spectra over several eV and 100 fs trajectories, considering couplings between all 80th lowest excited states, at low computational cost with a mixed semiempirical density functional theory setup.

<figure><img src="image_8.png" /></figure>

## 1. INTRODUCTION

The time-dependent study of photochemical properties in extended systems implying periodic boundary conditions, i.e., in solids, surfaces, or condensed-phase systems, requires a setup suitable for efficient nonadiabatic dynamics. Localization of the excitation manifold, ideally in combination with a hybrid quantum mechanics/molecular mechanics setup, can reduce computational timings and enable applications at a large scale.<sup>1,2</sup> In combination, trajectory surface hopping (TSH) represents a powerful tool to avoid the scanning of complete potential energy surfaces: Among the mixed quantum-classical approaches, TSH restricts the classical description of nuclear motion to sampling a single trajectory and thus enabling straightforward parallelization.<sup>3–9</sup> In general, internal consistency of trajectory surface hopping methods can be restored by adding decoherence corrections, then yielding satisfactory results in agreement with quantum wave packet computations.<sup>10</sup> The advantage of TSH becomes even more significant when modeling systems with periodic boundary conditions (PBC). In contrast to modeling distinct molecules, the number of relevant excited states is usually much larger for extended system sizes.<sup>11</sup> In these cases, TSH solely requires the computation of nuclear forces for the single currently

populated state, while Ehrenfest methods, e.g., refer to forces averaged across all states.<sup>8</sup>

Especially for molecular systems, the success of TSH can moreover be reasoned with further method developments, e.g., for the efficient computation of nonadiabatic couplings (NACs): Explicit analytical differentiation<sup>12–15</sup> is computationally demanding, and accurate integration of analytical couplings is still hampered due to the nevertheless required discretization in time.<sup>16,17</sup> Time derivative couplings<sup>17–19</sup> have therefore been suggested, relying on the numerical derivative of overlaps of wave functions of consecutive time steps. While the determinant derivative (DD) ansatz<sup>20,21</sup> refers to overlaps of Slater determinants, the orbital derivative (OD) ansatz restricts the contributions to overlap matrices of Kohn–Sham orbitals, reducing computation times by at least 2 orders of magnitude while showing similar accuracies as DD.<sup>22–25</sup> OD

Received: June 30, 2025

Revised: September 27, 2025

Accepted: October 1, 2025
requires tracking of the phase and energetic ordering of the underlying Kohn–Sham orbitals, and it has therefore been suggested to rely on unitary-transformed instead of canonical orbitals.²⁶ Furthermore, a modified OD approach applying an additional transformation step based on singular value decomposition has been developed to enhance numerical stability.²⁶,²⁷

Next to the mentioned ansätze for numerical time derivative couplings, local diabatization (LD) has been established as a standard approach for high stability in highly peaked nonadiabatic interactions, relying on overlap matrices computed within the DD or OD ansatz.²⁸,²⁹ Furthermore, the most cost-efficient alternative for computing NACs is given by phenomenological couplings, as defined, e.g., within the time-dependent Baeck-An (TDBA) approximation, relying solely on energy gaps and their second time derivatives.³⁰⁻³³ For ethylene and fulvene, it was shown that TDBA couplings provide a qualitatively correct description of the dynamics in comparison to analytical couplings, failing, however, in regions of strong nonadiabatic coupling and small velocities. More explicitly, adiabatic state populations and corresponding lifetimes are reproduced within the error bars of comparable results corresponding to exact couplings. The magnitudes of couplings as well as the corresponding order of coupling magnitudes for different pairs of states are, in general, well-reproduced with the trend to overestimate small-magnitude and underestimate large-magnitude couplings.

While the mentioned developments have significantly advanced TSH schemes in recent years, implementations and applications are often still restricted to molecular systems. Several trajectory surface hopping codes for periodic boundary conditions exist, including Pyxaid,³⁴ Libra,³⁵⁻³⁷ Hefeı,³⁸ CPMD,³⁹⁻⁴¹ ZAGREB,⁴² or GTSH,⁴³ enabling nonadiabatic molecular dynamic simulations through interfaces with state-of-the-art second-party electronic structure programs. However, many of these rely on the neglect of back reaction (NBRA) approximation or are restricted in terms of coupling schemes and modularity. NBRA still represents a common standard approximation, assuming that ground-state nuclear forces can be used to propagate the nuclei, approximating electronic excitation energies as simple Kohn–Sham energy gaps and time-derivative couplings as derivatives of those Kohn–Sham orbitals.⁴⁴,⁴⁵ The NBRA increases efficiency and is considered a reliable approximation in cases where electronic excitation is not associated with significant structural changes, which is often well justified for slow nonadiabatic processes, such as electron–hole recombination. The approximation is, however, inadequate in cases where coupling between nuclear and electronic structure needs to be taken into account, as has been reasoned, e.g., for an accurate modeling of charge transfer.⁴⁶ Furthermore, the mentioned interfaces for solid-state chemistry differ first of all in the underlying quantum-chemical electronic-structure program, providing thus only the corresponding characteristic features of these electronic-structure codes. Pyxaid and Hefeı-NAMD are currently interfaced with Quantum Espresso and VASP, respectively, relying thus on a global plane wave basis for the electronic-structure computation. CPMD comprises its own surface hopping module, which is equally restricted to plane waves. Regarding the underlying electronic-structure code, the closest comparison to the interface presented in this contribution is possible with the Libra, ZAGREB, and GTSH interfaces, which enable all computations using the CP2K program package. In

contrast to our contribution, the Libra-CP2K interface described in ref 47 employs the NBRA and calculates numerical time derivative couplings based on spin-adapted configurations, which represent Slater determinants for closed-shell structures. Combining CP2K with the Zagreb Surface Hopping NAMD code furthermore provides a setup for performing ΔSCF trajectory surface hopping, relying on Landau–Zener hopping probabilities.⁴² Finally, very recent work on the GTSH program package overcomes the NBRA, but is still restricted to OD couplings and assesses solely semiempirical GFN1-xTB-sTDA kernels for simulating photo-induced dynamics in carbon nanotubes.⁴³

To address these limitations, we present a flexible and efficient interface between the surface hopping package NEWTON-X⁴⁸⁻⁵⁰ and the electronic structure code CP2K,⁵¹ supporting periodic systems and multiple nonadiabatic coupling approaches. The electronic-structure setup relies on the mixed Gaussian and plane waves (GPW) framework⁵¹⁻⁵³ of CP2K (Section 2.1), which features both approximate hybrid density functional as well as semiempirical kernels for Sternheimer formulations within the Tamm-Dancoff approximation (TDA) of time-dependent density functional theory (TDDFT).⁵⁴⁻⁵⁸ NEWTON-X, in return, provides several choices for coupling constants (Section 2.2), including local diabatization,²⁸,²⁹ OD²³ and time-dependent Baeck-An couplings.³² The interface is assessed by comparing excited-state (ES) dynamics for molecular pyrazine with analogous simulations using the already established NEWTON-X interface with the molecular quantum chemistry program package TURBOMOLE⁵⁹,⁶⁰ (Section 3.1). Applicability and efficiency of the mixed density functional theory and semiempirical setup is furthermore assessed for crystalline pyrazine (Section 3.2).

## 2. METHODS

The underlying electronic-structure and surface hopping methods of the program packages CP2K⁵¹,⁵³⁻⁵⁸ and NEWTON-X³⁰,⁴⁸⁻⁵⁰ are well documented, and the following summary of computational models is thus restricted to summarizing the structure and versatility of the interface, as well as to provide explicit formulas for the computation of numerical orbital derivative couplings.

### 2.1. Approximate Hybrid Density Functional Accuracy for Adiabatic Wave Functions, Potential Energies, and Excited-State Forces.
Using the Born–Oppenheimer definition of adiabatic states, the total electronic wave function of the system $Ψ(\mathbf{R}(t))$ can be described as a linear expansion of such adiabatic eigenfunctions $ψ^M(\mathbf{R}(t))$ of all electronic states $M$ with corresponding expansion coefficients $d^M(t)$

$$|\Psi(\mathbf{R}(t))\rangle = \sum_{M} d^M(t)|\psi^M(\mathbf{R}(t))\rangle \quad (1)$$

depending on the nuclear coordinates $\mathbf{R}$ of time step $t$. The time propagation of coefficients $d^M(t)$ is hereby determined by the time-dependent Schrödinger equation

$$i \frac{d d^M(t)}{dt} = \sum_N d^N(t) [\delta_{MN} E_N(\mathbf{R}(t)) - i \tau_{MN}(t)] \quad (2)$$

referring to adiabatic potential energy surfaces $E_M(\mathbf{R}(t))$ and nonadiabatic coupling elements (NACs) $\tau_{MN}(t)$

$$\tau_{MN}(t) = \langle \psi^M(\mathbf{R}(t)) | \frac{\partial}{\partial t} \psi^N(\mathbf{R}(t)) \rangle \quad (3)$$
<figure><img src="image_9.png" /><figcaption>**Figure 1.** Interface of CP2K and NEWTON-X: NEWTON-X provides Cartesian coordinates of subsequent time steps $\mathbf{R}(t)$, $\mathbf{R}(t + \Delta t)$ to CP2K, which in return passes molecular orbital coefficients $C$, Cartesian normal modes $\omega$, excited-state eigenvalues $\Omega$ and amplitudes $\mathbf{X}$, oscillator strengths $f$, electronic-state forces F and, depending on the choice for coupling constants, the overlap matrix $\bar{S}$ over wave functions of two subsequent time steps. Next to the general infrastructure, including e.g., adiabatic dynamics, thermostats, and decoherence corrections, NEWTON-X enables choosing between time-dependent Baeck-An or orbital derivative couplings or to rely on local diabatization. CP2K offers the modular choice to combine various ground-state references, e.g., density functionals within the generalized gradient approximation (GGA), hybrid functionals implying the auxiliary density matrix method (ADMM) or semiempirical GFN1-xTB references with corresponding excited-state kernels within the simplified or conventional Tamm-Dancoff approximation, the former abbreviated as sTDA.</figcaption></figure>

The adiabatic excited-state wave functions $|\psi^M(\mathbf{R}(t))\rangle$ are not directly accessible within a time-dependent density functional theory framework, but can be approximated, implying a configuration interaction ansatz with singly excited Slater determinants $|\Phi_{ai\sigma}\rangle$

$$|\psi^M\rangle = \sum_{i\sigma} c_{i\sigma}^M |\Phi_{i\sigma}\rangle \quad (4)$$

The singly excited determinants $|\Phi_{i\sigma}\rangle$ describe an excitation from the occupied to the virtual molecular orbital space, indicated by indices $\{i, j, ...\}$ and $\{a, b, ...\}$, respectively. $\sigma$ denotes the electronic spin index. Within the Tamm-Dancoff approximation (TDA), the expansion coefficients $c^M$ are defined based on the excited-state TDA eigenvectors $\mathbf{X}^{M26,61}$

$$c_{i\sigma}^M = X_{i\sigma}^M \quad (5)$$

with the latter stemming from the Sternheimer formulation of the hermitian TDA eigenvalue problem

$$\mathbf{A}\mathbf{X}^M = \Omega_M \mathbf{S}\mathbf{X}^M \quad (6)$$

$$\sum_{kk} [F_{\mu\kappa\sigma}\delta_{ik} - F_{ik\sigma}S_{\mu\kappa}]X_{kk\sigma}^M + \sum_{\kappa\lambda} Q_{\mu\kappa}^T K_{\kappa\lambda\sigma} [\mathbf{X}^M] C_{\lambda i\sigma} \\ = \sum_{\kappa} \Omega_M S_{\mu\kappa} X_{\kappa i\sigma}^M \quad (7)$$

with corresponding excitation energies $\Omega_M = E_M - E_0$. $\mathbf{S}$ represents the standard overlap matrix over atomic orbitals, $\{\mu, \nu, ...\}$, and $\mathbf{Q}$ the projection operator onto the virtual orbital space inherent to the Sternheimer formulation. Within the latter, the virtual orbital space is omitted, thus the eigenvectors $X_{i\sigma}^M$ of eq 5 are transformed according to standard transformation rules from the molecular orbital to the atomic orbital basis with the MO coefficients $C$. Depending on the chosen ground-state reference and excited-state kernel, the explicit Kohn-Sham F and kernel K matrix elements comprise

different one- and two-electron contributions (see Figure 1); in the case of hybrid functionals, these would be one-electron, Coulomb, exact exchange, as well as functional-dependent exchange-correlation (XC) kernel contributions, with the computationally demanding exact exchange integrals being optionally approximated within the auxiliary density matrix method (ADMM).⁵⁷,⁶² Next to a conventional hybrid density functional kernel, semiempirical kernels based on the simplified Tamm-Dancoff approximation (sTDA) reduce computational timings, parametrizing Coulomb and exchange integral contributions.⁵⁷

We note that, within the TDA, hermiticity of the orbital-rotation matrix $\mathbf{A}$ ensures that the resulting CIS-type wave functions $\psi^M$ are orthogonal to each other and that eq 5 is an adequate choice for expansion coefficients $c^M$. Orthogonality is, however, in general not trivial and, when expanding the formalism to full TDDFT, mutually orthogonal wave functions need to be constructed by relying on the hermitian formulation of the full TDDFT eigenvalue problem.²⁶ Expansion coefficients $c^M$ are then given in terms of ground-state KS energy eigenvalues $\epsilon_{p\sigma}$ as⁶³

$$c_{i\sigma}^M = (\epsilon_{a\sigma} - \epsilon_{i\sigma})^{-1/2} (X_{i\sigma}^M + Y_{i\sigma}^M) \quad (8)$$

with the TDDFT deexcitation amplitudes $\mathbf{Y}^M$ being defined by solving the full TDDFT eigenvalue problem. Equation 8 thereby only holds for generalized gradient approximation (GGA) functionals, requiring further modifications for hybrid functionals, as explained in refs 63–65. Similar to eq 8, an alternative common choice for $c^M$ within TDA defines — in contrast to eq 5 — the coefficients in dependence of the excitation energy $\Omega_M$⁶⁶

$$c_{i\sigma}^M = \sqrt{\frac{\epsilon_{i\sigma} - \epsilon_{a\sigma}}{\Omega_M}} X_{i\sigma}^M \quad (9)$$

The following discussion will however always refer to eq 5.
Based on eq 7, ES nuclear gradients can be formulated by setting up a variational Lagrangian $L^M$ for the chosen excited state M, as outlined in ref 57. The force $F^M$, given as the negative derivative of the variational Lagrangian with respect to the nuclear coordinates R of nuclei $\alpha$

$$F_{\alpha}^{M}(t) = -\nabla_{\alpha} L^{M}(\mathbf{R}(t)) \qquad (10)$$

determines, according to Newton's first law of motion, the acceleration of each nuclei $\alpha$.

## 2.2. Efficient Approximations for Nonadiabatic Couplings.
Nonadiabatic coupling vectors (NACs) $\tau_{MN}(t)$ describe the interaction between electronic states M and N. Within Tully's fewest switches surface hopping,¹⁸ they are required for the computation of hopping probabilities. An analytical differentiation, defining NACs as the product of velocity $v$ times the analytical derivative of the adiabatic state $\psi^N$ with respect to the Cartesian nuclear coordinate

$$\begin{aligned} \tau_{MN}(t) &= \langle \psi^M(\mathbf{R}(t)) | \frac{\partial}{\partial t} \psi^N(\mathbf{R}(t)) \rangle \\ &= \mathbf{v} \cdot \langle \psi^M(\mathbf{R}(t)) | \nabla | \psi^N(\mathbf{R}(t)) \rangle \end{aligned} \qquad (11)$$

however, requires significant implementation and computational effort. Moreover, since analytical NACs are available at time step $t$ only, numerical differentiation, involving evaluation at two time steps $t$ and $t + \Delta t$, is often considered more robust. Several more cost-efficient approximations relying on the idea of numerical differentiation by Hammes-Schiffer and Tully¹⁸ have thus been established, approximating NACs with a 4-point formula according to

$$\tau_{MN}(t) \approx \frac{1}{4\Delta t} [3\bar{\tau}_{MN}(t, t - \Delta t) - 3\bar{\tau}_{NM}(t, t - \Delta t) \\ - \bar{\tau}_{MN}(t - \Delta t, t - 2\Delta t) + \bar{\tau}_{NM}(t - \Delta t, t - 2\Delta t)] \qquad (12)$$

$\bar{\tau}_{MN}(t, t - \Delta t)$ can be obtained either by straight-forwardly evaluating the resulting overlap matrices of Slater determinants (Determinant derivative (DD) ansatz²⁰,²¹,²⁴,²⁵,⁴⁹) or by first applying the Slater-Condon rules to reduce resulting contributions to overlap matrices of Kohn-Sham molecular orbitals (Orbital derivative (OD) approach²²,²³). More explicitly, DD overlap couplings are based on the wave function ansatz of eq 5 given as

$$\bar{\tau}_{MN}^{DD}(t, t - \Delta t) = \sum_{ia\sigma} X_{ia\sigma}^{M}(t) X_{jb\sigma}^{N}(t + \Delta t) \langle \Phi_{ia\sigma}(\mathbf{R}(t)) | \Phi_{jb\sigma}(\mathbf{R}(t + \Delta t)) \rangle \qquad (13)$$

implying the computation of overlap matrices over singly excited determinants $|\Phi_{ia\sigma}(\mathbf{R}(t))\rangle$. These overlap matrices can be evaluated by applying Löwdin's rule, requiring the computation of all pairs of possible determinants $|\Phi_{ia\sigma}(\mathbf{R}(t))\rangle$ and $|\Phi_{jb\sigma}(\mathbf{R}(t))\rangle$. In contrast, the OD ansatz is obtained by first applying the Slater-Condon rules and subsequently evaluating the resulting overlap expressions

$$\begin{aligned} \bar{\tau}_{MN}^{OD}(t, t - \Delta t) = & \sum_{ia\sigma} (X_{ia\sigma}^{M}(t))^T X_{ia\sigma}^{N}(t - \Delta t) \\ & + \sum_{ia\sigma} (X_{ia\sigma}^{M}(t))^T X_{ib\sigma}^{N}(t - \Delta t) \bar{S}_{ib\sigma}(t, t - \Delta t) \\ & - \sum_{ija\sigma} (X_{ia\sigma}^{M}(t))^T X_{aj\sigma}^{N}(t - \Delta t) \bar{S}_{j\sigma\sigma}(t, t - \Delta t) \end{aligned} \qquad (14)$$

The computational advantage results from the fact that $\bar{S}_{pq\sigma}(t, t - \Delta t)$ represents an overlap matrix over Kohn-Sham orbitals $|\phi_{p\sigma}(\mathbf{R}(t))\rangle$

$$\bar{S}_{pq\sigma}(t, t - \Delta t) = \langle \phi_{p\sigma}(\mathbf{R}(t)) | \phi_{q\sigma}(\mathbf{R}(t - \Delta t)) \rangle \qquad (15)$$

which is provided via the electronic-structure program. Further approximations have been suggested, with time-dependent Baeck-An couplings being a computationally efficient representative relying solely on electronic-state energies as input,³⁰

$$\begin{aligned} \tau_{MN}^{\text{TDBA}}(t, t - \Delta t) = & \left\{ \frac{\operatorname{sgn}(\Delta E_{MN})}{2} \sqrt{\frac{1}{\Delta E_{MN}}} \frac{d^2 \Delta E_{MN}}{dt^2} \right. & \text{if } \frac{1}{\Delta E_{MN}} \frac{d^2 \Delta E_{MN}}{dt^2} > 0 \\ & \left. 0 \right. & \text{if } \frac{1}{\Delta E_{MN}} \frac{d^2 \Delta E_{MN}}{dt^2} \le 0 \end{aligned} \qquad (16)$$

with $\Delta E_{MN} = E_M - E_N$. Last but not least, couplings become negligible when relying on diabatic wave functions, constructed via local diabatization⁵⁹ based on a transformation matrix that depends on the overlap matrix $\bar{S}_{pq\sigma}(t, t - \Delta t)$. Local diabatization interpolates between two time steps and thus avoids sampling of very peaked NACs, enabling the choice of larger time steps.

## 2.3. Interface Setup Including Initial Conditions and State and Phase Tracking.
A graphical overview of the technical architecture of the CP2K-NEWTON-X interface is given in Figure 1, detailing how inputs and outputs are exchanged and visualizing the modular option to choose different combinations of ground-state references, implying corresponding Kohn-Sham matrix elements $F$ and $M\Omega$ energies $\epsilon$, as well as excited-state kernels, with corresponding excitation energies $\Omega_M$ and excitation amplitudes $X^M$. As visualized in Figure 1, these quantities are parsed from CP2K to NEWTON-X and enable the computation of matching coupling constants. Data exchange thereby depends on whether the NEWTON-X classical series (CS) or NEWTON-X new series (NS) is employed, with the former being based on Perl scripts enabling data exchange via text files, while the latter relies on monolithic Fortran code to maximize computational performance and data management, now enabling data exchange through memory.⁵⁰ The interface is set up such that the complete general infrastructure of CP2K and NEWTON-X is made accessible, thus NEWTON-X is providing — next to the discussed features in the benchmark studies of this contribution — e.g., adiabatic dynamics, thermostats, and decoherence corrections. As outlined, choices for NACs comprise TDBA and OD couplings as well as local diabatization based on the latter. CP2K, in return, provides the modular choice to rely on various ground-state references, including density functionals within the generalized gradient approximation (GGA), hybrid functionals, also within the auxiliary density matrix method (ADMM), as well as semiempirical GFN1-xTB references.⁶⁷,⁶⁸ The thereon-based excited-state kernel can be chosen flexibly, thus enabling the combination of, e.g., GGA ground-state references with an ADMM-approximated hybrid functional excited-state kernel. Computational costs are furthermore reduced by 1 order of magnitude if relying on sTDA instead of the conventional TDA kernels.⁵⁷ Semiempirical approximations can thus be introduced for both the ground-state reference (within GFN1-xTB) as well as for the excited-state kernel (within sTDA). We
<figure><img src="image_10.png" /><figcaption>**Figure 2.** Left: Electronic absorption spectrum [in eV] of the lowest 5 excited states of pyrazine comparing implementations within the program packages TURBOMOLE and CP2K, using the B3LYP or PBE functional and an all-electron, molecular def2-TZVP or a MOOPT-TZVP basis set, the latter in combination with GTH pseudo potentials. Experimental reference data is taken from ref 70 Right: Exemplary trajectory for molecular pyrazine based on OD couplings, the B3LYP functional and a molecular def2-TZVP basis set. Energies are reported in atomic units (a.u.).</figcaption></figure>

also refer to refs 51,69 and 48–50 for further information on the CP2K and NEWTON-X program packages.

Furthermore, various options to set up an ensemble of initial states are automatically provided by NEWTON-X.⁵⁰ The nuclear ensemble approach, e.g., defines mass-scaled Cartesian coordinates and momenta based on a Wigner probability distribution $P_W$

$$P_W(\bar{\mathbf{R}}, \bar{\mathbf{P}}) = \prod_{I=1}^{N_F} \frac{\alpha_I}{\pi \hbar} \exp\left(-\frac{\alpha_I}{\hbar \omega_I} (\omega_I^2 \bar{R}_I^2 + \bar{P}_I^2)\right) \quad (17)$$

with $\alpha_I = \tanh(\frac{\hbar \omega_I}{2k_B T})$, $\bar{R}_I = \mu_I^{1/2} R_P$, $\bar{P}_I = \mu_I^{-1/2} P_P$, the reduced mass $\mu_I$ and vibrational frequencies $\omega_I$. For systems implying periodic boundary conditions lacking rotational degrees of freedom, the thermal population of vibrational states thereby refers to the reduced number of degrees of freedom $N_F = 3N_{at} - 3$. Starting from the so-generated manifold of initial structures, propagated geometries and velocities, $\mathbf{R}(t)$ and $\mathbf{v}(t)$, are obtained according to the classical equations of motion using the velocity Verlet algorithm

$$\mathbf{R}(t + \Delta t) = \mathbf{R}(t) + \mathbf{v}(t)\Delta t + \frac{1}{2}\mathbf{a}(t)\Delta t^2 \quad (18)$$

$$\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{1}{2}(\mathbf{a}(t) + \mathbf{a}(t + \Delta t))\Delta t \quad (19)$$

requiring the computation of acceleration $\mathbf{a}(t)$ as defined via eq 10 and thus the computation of both ES energies and ES nuclear gradients.

Furthermore, expansion coefficients $d_M(t)$ are propagated according to the time-dependent Schrödinger equation (eq 2) based on the computed NACs. The expansion coefficients $d_M(t)$ enable to propagate state occupations based on FSSH hopping probabilities $P$

$$P_{M \to N}(t) = \max \left[ 0, \frac{-2\Delta t}{|d_M(t)|^2} \tau_{NM}(t) \text{Re}[d_M(t) d_N^*(t)] \right] \quad (20)$$

i.e., to switch the state from $M$ to $N$ if the hopping probability fulfills $\sum_{K=1}^{N-1} P_{M \to K} < r_t \le \sum_{K=1}^{N} P_{M \to K}$ with $r_t$ representing a randomly chosen number drawn from a uniform distribution over the interval [0, 1]. To ensure correct time propagation of excited-state amplitudes and maintain continuity in state labeling, the interface enforces state tracking via overlap-based orthonormality checks between subsequent time steps

$$\sum_{ij\mu\nu\sigma} X_{\mu\nu\sigma}^{N}(t) S_{\mu\nu}(t) X_{\nu\sigma}^{M}(t - \Delta t) = \delta_{MN} \quad (21)$$

with the excited-state eigenvectors $X^M$ being again expressed within the AO-based formulation of the Sternheimer equations. Moreover, phase tracking is ensured by comparing the phase of the orthonormalized excited-state eigenvectors along trajectories.

### 3. RESULTS AND DISCUSSION

**3.1. Validation for Molecular Pyrazine.** Pyrazine represents an extensively studied molecule in photochemistry. The internal conversion from the second to the first bright excited state was investigated both experimentally with time-resolved photoelectron spectroscopy measurements, as well as theoretically using time-dependent density functional theory.⁴⁸,⁶⁴,⁷⁰⁻⁷⁴ The absorption spectrum of molecular pyrazine in the gas phase is thus characterized by two bands at energies of 3.8 and 4.8 eV, being classified as $n \to \pi^*$ transition (HOMO → LUMO) to the first singlet state of $^1B_{1u}$ character and a second $\pi \to \pi^*$ transition (HOMO-1 → LUMO) to the second excited singlet state of $^1B_{2u}$ symmetry. The third bright state corresponds to a $\pi \to \pi^*$ transition (HOMO-1 → LUMO+1) with $^1B_{3u}$ character. To validate the
implementation in comparison to the already interfaced molecular quantum chemistry code TURBOMOLE⁵⁹,⁶⁰ and to compare results with the extensive computational as well as experimental literature data, we computed the broadened absorption spectrum as well as the excited-state lifetimes for molecular pyrazine in the gas phase.

Computational setups for the molecule are chosen such that they enable a comparison with molecular quantum-chemistry codes, with all-electron computations being performed within the Gaussian and augmented plane wave framework (GAPW),⁵⁴,⁷⁵ relying on the B3LYP functional⁷⁶ in combination with the def2-TZVP basis set of Ahlrichs⁷⁷,⁷⁸ and all-electron pseudo potentials.⁷⁵ To further assess setups suitable for periodic boundary conditions, GAPW was replaced by the Gaussian and plane wave framework (GPW)⁷⁹ and all-electron basis sets by valence-only MOLOPT basis sets⁸⁰ in combination with corresponding pseudo potentials.⁸¹ Broadened absorption spectra were obtained using the nuclear ensemble approach, sampling over 500 initial geometries and velocities generated according to a Wigner probability distribution.⁸² For the simulated cross section, a Lorentzian line shape with a broadening of 0.05 eV was used. Data files relevant for the presented computational results are provided in ref 83 and on Materials Cloud.

The broadened absorption spectrum is shown in Figure 2; excitation energies of the lowest states corresponding to excitations from HOMO→LUMO (B₁u), HOMO→LUMO+1 (A_u), and HOMO−1→LUMO (B₂u) are given in Table 1. For

<table><caption>**Table 1. Excitation Energies [in eV] for the First Three Excited States of ¹B₁u (HOMO→LUMO), ¹A_u (HOMO→LUMO+1), and ¹B₂u (HOMO−1→LUMO) of Molecular Pyrazine at Various Levels of Theory**</caption><thead><tr><th></th><th>setup</th><th>B<sub>1u</sub> S<sub>1</sub></th><th>A<sub>u</sub> S<sub>2</sub></th><th>B<sub>2u</sub> S<sub>3</sub> or S<sub>4</sub></th></tr></thead><tbody><tr><td>Exp.<sup>64</sup></td><td></td><td>3.8</td><td>—</td><td>5.69</td></tr><tr><td>TURBOMOLE</td><td>B3LYP/TZVP</td><td>4.00</td><td>4.62</td><td>S<sub>3</sub>: 5.61</td></tr><tr><td>CP2K</td><td>GAPW/B3LYP/TZVP</td><td>4.00</td><td>4.62</td><td>S<sub>3</sub>: 5.59</td></tr><tr><td></td><td>GPW/B3LYP/MOL-TZ</td><td>4.07</td><td>4.76</td><td>S<sub>3</sub>: 5.60</td></tr><tr><td></td><td>GPW/PBE/MOL-TZ</td><td>3.60</td><td>4.05</td><td>S<sub>4</sub>: 5.42</td></tr><tr><td></td><td>sTDA/PBE/MOL-DZ</td><td>3.27</td><td>3.95</td><td>S<sub>4</sub>: 5.32</td></tr></tbody></table>

the sake of convenience, these states will be referred to as S₁, S₂, and S₃, with the index referring to the energetic ordering of theoretical results. Within GAPW and relying on the B3LYP method, the CP2K results reproduce the TURBOMOLE reference data, validating the implementation of the interface. Switching from molecular basis sets to MOLOPT basis sets, optimized for PBC, introduces only small deviations of max. 0.07 eV. Relying, however, on the GGA functional PBE⁸⁴ or switching to sTDA kernels reorders the excited states, with an excitation from HOMO−2→LUMO being lower in energy than the excitation to the B₂u state. The latter now corresponds to the fourth excited state with an excitation energy which is, by fortuitous error cancellation, closer to the experimental value, but differs from the B3LYP reference by 0.2−0.4 eV.

Excited-state relaxation is therefore started from state S₄ for PBE and sTDA@PBE kernels. Furthermore, the first B₁u and the dark A_u state are shifted to lower excitation energies of 3.60/3.27 and 4.05/3.95 eV for PBE and sTDA@PBE kernels, respectively. Further setups were tested, and corresponding results are given in the Supporting Information (Figure S2).

The fact that GGA functionals reverse the energetic ordering of excited states in comparison to hybrid functionals underlines the importance of accurate XC functionals to describe excited states.⁸⁵ It might thus be beneficial to test optimally tuned range-separated hybrid functionals, which outperform hybrids and were shown to have an accuracy for energy gaps, excited-state characters, and NACs comparable to semiempirically parametrized CAM-B3LYP or EOM-CCSD.⁸⁶

To simulate excited-state relaxation, we performed non-adiabatic dynamics propagation of at least 200 trajectories for 200 fs, comparing time-dependent excited-state populations and lifetimes. An exemplary trajectory is visualized in Figure 2, highlighting that, starting from the second bright state ¹B₂u relaxation to the lower dark state ¹A_u is observed already within the first 10 fs, subsequently followed by a transition to the energetically lowest bright ¹B₁u state. The latter is well separated from the ground state, so that no relaxation to the ground state is observed. Corresponding adiabatic populations are visualized in Figure 3, comparing the various choices for coupling constants (OD, TDBA, and LD with OD overlap matrix) that can be combined flexibly with the various quantum-mechanical setups (based on GAPW and GPW in combination with conventional TDA kernels or the simplified Tamm-Dancoff approximation (sTDA)). Additional plots further validating, e.g., the agreement between TURBOMOLE and CP2K or adding further setups with varying coupling constant - quantum mechanical method combinations are given in the Supporting Information (Figure S3). Populations plotted in Figure 3a, comparing CP2K and TURBOMOLE results for Baeck-An couplings, again, validate the interface with respect to the existing setups.

Furthermore, comparing the different coupling schemes, state populations agree well. However, even though the fitted lifetime constants of approximately 10 fs agree well among the different setups, they are too small in comparison to the literature values of 20−22 fs.⁶⁴,⁸⁷−⁸⁹ Comparing the computational setup of the TURBOMOLE computations with ref 64 suggests that the main source of error is to be seen in the underlying coupling scheme, with DD couplings thus outperforming TDBA or OD couplings. Switching from conventional TDA to sTDA kernels (see Figure 3c,d) includes, due to starting from S₄, subsequent population of S₃ and S₂, with similar population profiles for both OD and TDBA couplings, as well as local diabatization. Relaxation is, however, even more overestimated with lifetimes of approximately 5-12 fs. For all density-functional theory and semiempirical setups, relaxation to the ground state is not observed within the first 200 fs.⁹⁰

Coupling magnitudes have been compared for OD couplings for CP2K and TURBOMOLE implementations, plotting the absolute coupling magnitude along a trajectory of 100 fs for the coupling of the two lowest excited states S₁ and S₂ as well as for the subsequent higher-lying excited states S₂ and S₃ for two different trajectories. The match in coupling magnitudes again validates the agreement of both implementations. A comparison of coupling magnitudes between different coupling setups, e.g., for OD vs TDBA couplings, is moreover given in Figure 4, providing an overall qualitative agreement, with the time and width of the couplings being well-reproduced. Absolute magnitudes are, however, underestimated by TDBA. The latter could be due to the chosen default setup for the TDBA parameters, with TDBA couplings being set to zero if the absolute energy gap of the involved states, $δε = |ΔE_{MN}(t)|$, was larger than 2 eV. Benchmark results of ref 30 suggest that an
<figure><img src="image_11.png" /><figcaption>**Figure 3.** Time-dependent excited-state adiabatic populations as obtained by averaging over 200 nonadiabatic trajectories for pyrazine, validating the CP2K implementation with respect to TURBOMOLE, and comparing 1. the different coupling setups — (a) and (c) TDBA, (b) OD, (d) LD with OD overlap — as well as 2. the different quantum-mechanical kernels — (a) and (b) conventional GAPW, (c) and (d) simplified Tamm-Dancoff approximation (sTDA). All computations were performed using the B3LYP functional and a def2-TZVP basis set.</figcaption></figure>

improvement in coupling intensities could be achieved by optimizing those parameters. It thus has to be stressed that Figure 4 can not be referred to when aiming for a general discussion on the reliability of TDBA for molecular systems. Restricting the viewpoint to setting a nonoptimized reference for the comparison with results for crystalline pyrazine, implying periodic boundary conditions, this nonoptimized choice of parameters for TDBA couplings is, however, sufficient. Further benchmarking to determine an optimal

parameter choice is not the focus of the current contribution and is therefore not discussed further.

**3.2. Extension to Crystalline Pyrazine.** Pure crystalline pyrazine does not show any fluorescence.⁹¹ Experimental spectroscopic data on crystalline pyrazine are restricted to phosphorescence,⁹¹,⁹² competing with quenching via triplet-triplet annihilation and showing lifetimes of the lowest triplet state in the range of 10⁻² s. However, pyrazine-based dyes with modified functional groups are broadly used to tune fluorescence in the solid state,⁹³⁻⁹⁶ with aggregation-induced
<figure><img src="image_12.png" /><figcaption>**Figure 4.** Comparison of coupling magnitudes $|σ_{ij}|$ [in a.u.] for OD and TDBA couplings along two exemplary trajectories of 100 fs for molecular pyrazine, visualizing couplings between the lowest $B_{1u}$ and the dark $A_u$ state (OD in orange, TDBA in blue) as well as the couplings between the second $B_{2u}$ and the dark $A_u$ state (OD in light green, TDBA in dark green).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 5.** Left: Electronic absorption spectrum of the lowest 80 excited states of crystalline pyrazine, using the B3LYP or PBE functional or sTDA@PBE kernels, a MOLOPT-TZVP or MOLOPT-DZVP basis set and GTH pseudo potentials. Experimental reference data⁹⁸ discriminates between the absorption spectrum obtained with light incident perpendicular to the (101) crystal face (Exp $⊥$), spectra obtained from projection onto the (100) crystal face (Exp $∥$) as well as the solution spectrum in hexane (Exp Solution). Right: Exemplary trajectory for crystalline pyrazine based on orbital derivative couplings, an sTDA kernel, the PBE functional, and a MOLOPT-DZVP basis set. Excited states of higher energy are so dense that single trajectories cannot be distinguished. To nevertheless associate bands with a range of excited states, color coding is adjusted, always highlighting 20 excited states with the same color.</figcaption></figure>

emission representing, e.g., a versatile tool to modify photophysical properties.⁹⁷ Solid-state absorption and fluorescence properties were hereby found to be more sensitive to (electron-donating) abilities of functional groups due to crystallization. An extensive study on tuning triplet-triplet annihilation or exploiting aggregation-induced emission for

different materials classes, such as functionalized pyrazine dyes, is beyond the scope of this contribution. For the current study, we limit the discussion to the broadened absorption spectrum of crystalline pyrazine in comparison to experimental reference data⁹⁸ and demonstrate the usability and efficiency regarding computational timings of the outlined tool set.
<figure><img src="image_14.png" /><figcaption>**Figure 6.** Time-dependent adiabatic populations as obtained by averaging over 91 nonadiabatic trajectories for crystalline pyrazine, relying on (a) OD, (b) phenomenological time-dependent Baeck-An couplings or (c) local diabatization, sTDA, the PBE functional and a MOLOPT-DZVP basis set. TDBA leads to the artifact that the ground state is populated, with the corresponding population thus being additionally given in subfigure (b).</figcaption></figure>

Computational setups for the crystal are now consistently relying on the Gaussian and plane wave (GPW) method, using MOLOPT basis sets⁸⁰ and Goedecker-Teter-Hutter pseudopotentials.⁸¹,⁹⁹ 80 excited states were considered for the simulations in total. While different basis set sizes and functionals were tested for the computation of the broadened absorption spectrum, the setup was restricted, for the sake of convenience, for the surface hopping computations to sTDA kernels and relying on PBE⁸⁴ references generated with MOLOPT-DZVP basis sets. In this case, 91 surface hopping trajectories of 100 fs were computed. While this sample size is

smaller than the 200 trajectories run for molecular pyrazine, simulations were checked to be converged with respect to the number of trajectories. To avoid artifacts due to the high density of states,¹⁰⁰ the nuclear time step was reduced to 0.1 fs in comparison to molecular pyrazine. All pairwise couplings between all 80 excited states were computed at each time step rather than restricting to adjacent-state interactions. The crystal structure is visualized in the Supporting Information (Figure S1) with the unit cell comprising 40 atoms.

Figure 5 visualizes on the left-hand side the broadened absorption spectrum for the various computational setups as
<figure><img src="image_15.png" /><figcaption>**Figure 7.** Comparison of coupling magnitudes $|σ_{ij}|$ [in a.u.] for OD and TDBA couplings along a typical trajectory of 100 fs for crystalline pyrazine, visualizing couplings between excited states S₁ and S₂ (left-hand side; OD in light green, TDBA in dark green) as well as the couplings between excited states S₂ and S₃ (right-hand side; OD in orange, TDBA in blue).</figcaption></figure>

well as experimental measurements. Plotted experimental data were taken from ref 98 and discriminate between the absorption spectrum obtained with light incident perpendicular to the (101) crystal face (denoted as Exp $⊥$), spectra obtained from projection onto the (100) crystal face (denoted as Exp $∥$) as well as the solution spectrum in hexane. The experimental spectrum thereby exhibits two distinct broad peaks: a smaller one from 3.7 to 4.2 eV and a larger one from 4.2 to 5.2 eV, peaking at 4.7 eV. The computational results for the 80 lowest states cover an energy range of 2 to 7 eV, with the relatively broad absorption peak being centered for the PBE functional between 4.7 to 4.9 eV, only slightly depending on the chosen basis set. The center of the analogous B3LYP band is shifted to 5.5 eV and moreover shows a shoulder at 4.9 eV. Switching from conventional to sTDA kernels has a negligible effect on the overall spectrum with similar band widths and heights.

When comparing energies of the highest 10 excited states (71st to 80th), the energy differences nevertheless amount to 0.5 eV. Since the first benchmark results for molecular pyrazine indicated that semiempirical GFN1-xTB ground-state references need to be tuned, i.e., by shifting the virtual orbital space (see Figure S2 of the Supporting Information), an extension of the benchmark to further modify the underlying ground-state reference is beyond the scope of the current investigation. It should, however, be stressed that optimized GFN1-xTB references could further increase the computational efficiency by an order of magnitude and that the suitability of such setups was already tested and approved.³⁵,³⁶,⁴³ The band peak for sTDA kernel and PBE references corresponds to the 60th or 61st excited state, with both states being thus used to generate starting velocities and geometries, and the manifold of trajectories to obtain averaged adiabatic populations thus relies on the inputs obtained from these two states.

An exemplary trajectory for OD couplings is visualized in Figure 5 on the right-hand side, highlighting the increased

density of states as well as the fast relaxation, consecutively over the intermediate excited states, to the first excited state S₁ within the first 40 fs. Similar plots are obtained for LD couplings.

Resulting adiabatic populations, obtained by averaging over 91 trajectories, confirm this finding (Figure 6): All trajectories exhibit ultrafast internal conversion toward the first excited state S₁ with most population being transferred within 40–50 fs. For OD couplings, the high-lying states are quickly depopulated and, in return, occupation of increasingly lower excited states sets in, with a maximum population of the 20th excited state being found after ~8 fs and of the 10th excited state after ~25 fs, respectively. The population of the three lowest excited states increases after 20 fs, with the third and second excited states reaching a maximum at 40 fs and then slowly decreasing to below 5%. While the first excited state reaches a maximum population of 95%, no relaxation to the ground state is observed within the total simulation time of 100 fs. LD based on the OD overlap matrix yields a similar picture: After 100 fs, the first excited state reaches an average population of 90%, with the relaxation proceeding over lower lying states reflected in the subsequent maxima in respective populations.

Regarding analogous simulations relying on time-dependent Baeck-An couplings, it needs to be stressed that, as discussed within ref 30, the model is formally only valid to describe the coupling of two states. Computations were thus performed based on the assumption that two-state coupling is adequate also for crystalline pyrazine. Due to the much higher density of states, the TDBA input parameter $δε$ had to be adjusted accordingly to yield similar relaxation profiles within the first femtoseconds. Best agreement was found when setting TDBA couplings to zero if the absolute energy gap of the involved states, $δε$, was larger than 0.5 eV. As documented in Figure S5 of the Supporting Information, other variations of $δε$ result in increasing or decreasing relaxation times; further modification
of input parameters, to, e.g., guarantee smoothness in the variation of the couplings, did not improve the results. As visualized in graph b of Figure 6, relaxation is then comparable to OD couplings for the first 40 fs. After 46 fs, however, the electronic-state population of the ground-state increases above 0.1%, with steadily decreasing population of the three lowest excited states.

While the increasing population of the ground state flaws the performance of TDBA couplings, a comparison of coupling magnitudes for OD and TDBA couplings, as given in Figure 7, nevertheless reveals the overall quantitative agreement of both ansätze when restricting the assessment to internal conversion within the manifold of excited states (see also additional data of Figure S4 in the Supporting Information). Couplings between excited states S₁ and S₂ are consistently smaller than 0.02 a.u for both coupling models. In contrast, the interaction between excited states S₂ and S₃ is for both OD and TDBA systematically larger, even though TDBA underestimates corresponding OD coupling magnitudes.

In total, one can thus conclude that OD and LD yield similar population patterns, with S₁ reaching a population of ≈90–95%. TDBA couplings, in contrast, exhibit spurious decay to the ground state, even though a comparison of coupling magnitudes for relaxation within the excited-state manifold provides qualitative agreement with OD couplings.

The CP2K-NEWTON-X interface is set up with the motivation to provide sufficient computational efficiency for large-scale computations on extended systems. Based on the outlined preliminary assessment of different setups for molecular and crystalline pyrazine, computational timings representing the thereon based main result of this manuscript are summarized in Table 2 for crystalline pyrazine, referring

<table><caption>Table 2. Timings for Performing a NAMD Step for Crystalline Pyrazine Considering the Lowest 80 Excited States<sup>a</sup></caption><thead><tr><th rowspan="2">timing</th><th>TDA</th><th>sTDA</th><th>sTDA</th></tr><tr><th>PBE/MOL-DZ</th><th>PBE/MOL-DZ</th><th>GFN1-xTB</th></tr></thead><tbody><tr><td>CP2K</td><td>3 min 55 s</td><td>30 s</td><td>10 s</td></tr><tr><td>OD overlap</td><td>1 min 2 s</td><td>1 min 2 s</td><td>11 s</td></tr><tr><td>total OD</td><td>5 min 25 s</td><td>1 min 57 s</td><td>45 s</td></tr><tr><td>total TDBA</td><td>5 min 4 s</td><td>1 min 43 s</td><td>38 s</td></tr><tr><td>total LD</td><td>5 min 35 s</td><td>1 min 35 s</td><td>23 s</td></tr></tbody></table>

again to a unit cell comprising 40 atoms corresponding to 416 basis functions for a MOLOPT-DZ basis. Computations were performed on an AMD EPYC GENOA 9534 processor utilizing a total of 64 cores. For the electronic structure computation within CP2K computing 80 excitation energies as well as the gradient for the 61st excited state, the total computation time amounts to 30 s for MOLOPT-DZ basis sets, PBE references, and sTDA kernels, assuming tight convergence criteria of 10⁻⁹ a.u., as well as grid cutoffs of 800 Ry. Computing corresponding nonadiabatic coupling constants within the orbital derivative ansatz adds another minute to the total computation time of 1 min 57 s, demonstrating that the method is tractable even for crystalline systems. Relying on TDBA couplings or local diabatization yields comparable computation times of 1 min 43 s and 1 min 35 s, respectively. Computational costs increase by a factor of 5

when switching to conventional TDA kernels using GGA functionals and are reduced by a factor of 2 for semiempirical GFN1-xTB references.

Further computational timings for the underlying electronic structure computation were already reported in ref 57. It highlights that sTDA in general enables an order of magnitude in speedup compared to ADMM-approximated hybrid functional kernels. Thus, it is most beneficial when aiming to compute broadband spectra, which implies a large number of excited states. Relying on sTDA kernels based on GGA functionals, even larger unit cell sizes of hundreds of atoms should thus be accessible at a moderate computational cost, with timings in the order of minutes to compute both broadband electronic excitation spectra as well as the corresponding excited-state nuclear gradient.

## 4. CONCLUSIONS

With the outlined variety of time-dependent density functional theory and semiempirical approximations for electronic-structure methods and the range of time-dependent Baeck-An to numerical couplings for the trajectory surface hopping, the current contribution establishes the interface of the electronic-structure program package CP2K<sup>51,53</sup> and the surface-hopping code NEWTON-X<sup>48–50</sup> as a versatile tool set for nonadiabatic molecular dynamics on extended materials, implying periodic boundary conditions.

Assessed features from the electronic-structure program package comprise conventional TDDFT based on GGA as well as hybrid functionals, the latter accelerated with ADMM. Semiempirical sTDA complements the electronic-structure setup, which can be used relying on conventional DFT references or on semiempirical GFN1-xTB references, with the latter, however, requiring further adjustments and benchmarking. Both excited-state kernels can be combined either with GPW or GAPW. NEWTON-X modularly provides phenomenological time-dependent Baeck-An, numerical orbital time-derivative couplings, or enables the use of local diabatization based on the orbital derivative overlap matrix.

Validation of the interface was demonstrated for molecular pyrazine, highlighting the close comparison to the molecular electronic structure code TURBOMOLE<sup>59,60</sup> and demonstrating the adequacy of the mentioned approximations for coupling constants to describe excited-state relaxation from the second to the first excited state. The obtained results reflect the error range of the density functional theory setup, with excitation bands being off by 0.2 to 0.8 eV for the investigated lowest excited states, respectively. Lifetimes of 5–12 fs underestimate the experimental reference value of 22 fs.

Analogous computations for crystalline pyrazine underscore the agreement between lifetimes obtained for OD couplings and corresponding results relying on local diabatization, demonstrating the efficiency of the interface. This enables the performance of a NAMD step with a semiempirical TDDFT kernel, a DFT ground-state reference, and a double-ζ basis set within minutes. While both local diabatization as well as OD and TDBA couplings yield similar results for molecular pyrazine, care must be taken when using TDBA couplings for crystalline pyrazine, implying periodic boundary conditions and a comparably high density of states. In this case, the two-state nature of the TDBA model leads to population leakage into the ground state. Future work will be focused on extending the tool set to enable a more extensive and accurate

<footnote><sup>a</sup>Computations were performed on an AMD EPYC GENOA 9534 processor using 64 cores.</footnote>
description of excited-state phenomena, including e.g.,
intersystem crossing.

■ ASSOCIATED CONTENT

Supporting Information

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acs.jctc.5c01082.

Structural information on pyrazine crystal; visualization
of structural data broadened absorption spectra for
molecular pyrazine; adiabatic populations for molecular
pyrazine; coupling magnitudes for crystalline pyrazine;
adiabatic populations for TDBA couplings for crystalline
pyrazine (PDF)

Additional information (ZIP)

■ AUTHOR INFORMATION

Corresponding Author

**Anna Hehn** – Christian-Albrechts-University Kiel, 24118 Kiel, Germany; Kiel Nano, Surface and Interface Science, 24118 Kiel, Germany; ● orcid.org/0000-0002-0498-740X; Email: hehn@pctc.uni-kiel.de

Authors

Jan-Robert Vogt – Christian-Albrechts-University Kiel, 24118 Kiel, Germany

Michael Schulz – Christian-Albrechts-University Kiel, 24118 Kiel, Germany

Rafael Souza Mattos – Aix Marseille University, CNRS, ICR, 13397 Marseille, France; ● orcid.org/0000-0003-0215-7100

Mario Barbatti – Aix Marseille University, CNRS, ICR, 13397 Marseille, France; Institut Universitaire de France, 75231 Paris, France; ● orcid.org/0000-0001-9336-6607

Maurizio Persico – Dipartimento di Chimica e Chimica Industriale, S6124 Pisa, Italy; ● orcid.org/0000-0001-6861-9289

Giovanni Granucci – Dipartimento di Chimica e Chimica Industriale, S6124 Pisa, Italy; ● orcid.org/0000-0002-4753-6318

Jürg Hutter – Department of Chemistry, University of Zurich, 8057 Zurich, Switzerland

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jctc.5c01082

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

R.S.M. and M.B. received support from the French govern-
ment under the France 2030 investment plan as part of the
Initiative d'Excellence d'Aix-Marseille Université (AMIDEX
AMX-22-IN1-48, MLCHEM) and from the European
Research Council (ERC) Advanced Grant SubNano (grant
agreement 832237).

■ REFERENCES

(1) Groenhof, G. *Biomolecular Simulations: Methods and Protocols*; Monticelli, L.; Salonen, E., Eds.; Humana Press: Totowa, NJ, 2013; pp 43–66.

(2) Avagliano, D.; Bonfanti, M.; Garavelli, M.; González, L. QM/MM Nonadiabatic Dynamics: the SHARC/COBRAMM Approach. *J. Chem. Theory Comput.* **2021**, *17*, 4639–4647.

(3) Wang, L.; Qiu, J.; Bai, X.; Xu, J. Surface hopping methods for nonadiabatic dynamics in extended systems. *WIREs Comput. Mol. Sci.* **2020**, *10*, No. e1435.

(4) Lei, Y.; Ma, H.; Vasquez, L. Ab initio nonadiabatic dynamics of semiconductor materials via surface hopping method†. *Chin. J. Chem. Phys.* **2022**, *35*, 16–37.

(5) Xie, B.-B.; Jia, P.-K.; Wang, K.-X.; Chen, W.-K.; Liu, X.-Y.; Cui, G. Generalized Ab Initio Nonadiabatic Dynamics Simulation Methods from Molecular to Extended Systems. *J. Phys. Chem. A* **2022**, *126*, 1789–1804.

(6) Akimov, A. V.; Neukirch, A. J.; Prezhdo, O. V. Theoretical Insights into Photoinduced Charge Transfer and Catalysis at Oxide Interfaces. *Chem. Rev.* **2013**, *113*, 4496–4565.

(7) Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum–Classical Dynamics. *Chem. Rev.* **2018**, *118*, 7026–7068.

(8) Nelson, T. R.; White, A. J.; Bjorgaard, J. A.; Sifain, A. E.; Zhang, Y.; Nebgen, B.; Fernandez-Alberti, S.; Mozyrsky, D.; Roitberg, A. E.; Tretiak, S. Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials. *Chem. Rev.* **2020**, *120*, 2215–2287.

(9) Smith, B.; Akimov, A. V. Modeling nonadiabatic dynamics in condensed matter materials: some recent advances and applications. *J. Phys.:Condens. Matter* **2020**, *32*, No. 073001.

(10) Granucci, G.; Persico, M. Critical appraisal of the fewest
switches algorithm for surface hopping. *J. Chem. Phys.* **2007**, *126*,
No. 134114.

(11) Liu, D.; Wang, B.; Vasenko, A. S.; Prezhdo, O. V. Decoherence ensures convergence of non-adiabatic molecular dynamics with number of states. *J. Chem. Phys.* **2024**, *161*, No. 064104.

(12) Send, R.; Furche, F. First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance. *J. Chem. Phys.* **2010**, *132*, No. 044107.

(13) Ou, Q.; Bellchambers, G. D.; Furche, F.; Subotnik, J. E. First-order derivative couplings between excited states from adiabatic TDDFT response theory. *J. Chem. Phys.* **2015**, *142*, No. 064114.

(14) Tapavicza, E.; Tavernelli, I.; Rothlisberger, U. Trajectory Surface Hopping within Linear Response Time-Dependent Density-Functional Theory. *Phys. Rev. Lett.* **2007**, *98*, No. 023001.

(15) Tavernelli, I.; Tapavicza, E.; Rothlisberger, U. Nonadiabatic coupling vectors within linear response time-dependent density functional theory. *J. Chem. Phys.* **2009**, *130*, No. 124107.

(16) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Artifacts due to trivial unavoided crossings in the modeling of photoinduced energy transfer dynamics in extended conjugated molecules. *Chem. Phys. Lett.* **2013**, *590*, 208–213.

(17) Meek, G. A.; Levine, B. G. Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations. *J. Phys. Chem. Lett.* **2014**, *5*, 2351–2356.

(18) Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. *J. Chem. Phys.* **1994**, *101*, 4657–4667.

(19) Shu, Y.; Zhang, L.; Sun, S.; Truhlar, D. G. Time-Derivative Couplings for Self-Consistent Electronically Nonadiabatic Dynamics. *J. Chem. Theory Comput.* **2020**, *16*, 4098–4106.

(20) Pittner, J.; Lischka, H.; Barbatti, M. Optimization of mixed quantum-classical dynamics: Time-derivative coupling terms and selected couplings. *Chem. Phys.* **2009**, *356*, 147–152.

(21) Plasser, F.; Ruckenbauer, M.; Mai, S.; Oppel, M.; Marquetand, P.; González, L. Efficient and Flexible Computation of Many-Electron Wave Function Overlaps. *J. Chem. Theory Comput.* **2016**, *12*, 1207–1219.

(22) Ryabinkin, I. G.; Nagesh, J.; Izmaylov, A. F. Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum-Classical Methods. *J. Phys. Chem. Lett.* **2015**, *6*, 4200–4203.

(23) Stojanović, L.; Bai, S.; Nagesh, J.; Izmaylov, A. F.; Crespo-Otero, R.; Lischka, H.; Barbatti, M. New Insights into the State Trapping of UV-Excited Thymine. *Molecules* **2016**, *21*, 1603.
(24) Lee, S.; Kim, E.; Lee, S.; Choi, C. H. Fast Overlap Evaluations for Nonadiabatic Molecular Dynamics Simulations: Applications to SF-TDDFT and TDDFT. *J. Chem. Theory Comput.* **2019**, *15*, 882–891.

(25) Lee, S.; Horbatenko, Y.; Filatov, M.; Choi, C. H. Fast and Accurate Computation of Nonadiabatic Coupling Matrix Elements Using the Truncated Leibniz Formula and Mixed-Reference Spin-Flip Time-Dependent Density Functional Theory. *J. Phys. Chem. Lett.* **2021**, *12*, 4722–4728.

(26) Uratani, H.; Morioka, T.; Yoshikawa, T.; Nakai, H. Fast Nonadiabatic Molecular Dynamics via Spin-Flip Time-Dependent Density-Functional Tight-Binding Approach: Application to Non-radiative Relaxation of Tetraphenylethylene with Locked Aromatic Rings. *J. Chem. Theory Comput.* **2020**, *16*, 7299–7313.

(27) Uratani, H.; Yoshikawa, T.; Nakai, H. Trajectory Surface Hopping Approach to Condensed-Phase Nonradiative Relaxation Dynamics Using Divide-and-Conquer Spin-Flip Time-Dependent Density-Functional Tight Binding. *J. Chem. Theory Comput.* **2021**, *17*, 1290–1300.

(28) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. *J. Chem. Phys.* **2001**, *114*, 10608–10615.

(29) Plasser, F.; Granucci, G.; Pittner, J.; Barbatti, M.; Persico, M.; Lischka, H. Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer. *J. Chem. Phys.* **2012**, *137*, No. 22A514.

(30) do Casal, M. T.; Toldo, J. M.; Pinheiro, M.; et al. Fewest switches surface hopping with Baeck-An couplings. *Open Res. Europe* **2022**, 1, 49.

(31) Baeck, K. K.; An, H. Practical approximation of the non-adiabatic coupling terms for same-symmetry interstate crossings by using adiabatic potential energies only. *J. Chem. Phys.* **2017**, *146*, No. 064107.

(32) An, H.; Baeck, K. K. Practical and reliable approximation of nonadiabatic coupling terms between triplet electronic states using only adiabatic potential energies. *Chem. Phys. Lett.* **2018**, *696*, 100–105.

(33) Shu, Y.; Zhang, L.; Chen, X.; Sun, S.; Huang, Y.; Truhlar, D. G. Nonadiabatic Dynamics Algorithms with Only Potential Energies and Gradients: Curvature-Driven Coherent Switching with Decay of Mixing and Curvature-Driven Trajectory Surface Hopping. *J. Chem. Theory Comput.* **2022**, *18*, 1320–1328.

(34) Akimov, A. V.; Prezhdo, O. V. Advanced Capabilities of the PYXAID Program: Integration Schemes, Decoherence Effects, Multiexcitonic States, and Field-Matter Interaction. *J. Chem. Theory Comput.* **2014**, *10*, 789–804.

(35) Shakiba, M.; Smith, B.; Li, W.; Dutra, M.; Jain, A.; Sun, X.; Garashchuk, S.; Akimov, A. Libra: A modular software library for quantum nonadiabatic dynamics. *Software Impacts* **2022**, *14*, No. 100445.

(36) Smith, B.; Shakiba, M.; Akimov, A. V. Crystal Symmetry and Static Electron Correlation Greatly Accelerate Nonradiative Dynamics in Lead Halide Perovskites. *J. Phys. Chem. Lett.* **2021**, *12*, 2444–2453.

(37) Smith, B.; Shakiba, M.; Akimov, A. V. Nonadiabatic Dynamics in Si and CdSe Nanoclusters: Many-Body vs Single-Particle Treatment of Excited States. *J. Chem. Theory Comput.* **2021**, *17*, 678–693.

(38) Zheng, Q.; Chu, W.; Zhao, C.; Zhang, L.; Guo, H.; Wang, Y.; Jiang, X.; Zhao, J. Ab initio nonadiabatic molecular dynamics investigations on the excited carriers in condensed matter systems. *WIREs Comput. Mol. Sci.* **2019**, 9, No. e1411.

(39) Franco de Carvalho, F.; Tavernelli, I. Nonadiabatic dynamics with intersystem crossings: A time-dependent density functional theory implementation. *J. Chem. Phys.* **2015**, *143*, No. 224105.

(40) Hutter, J. Excited state nuclear forces from the Tamm-Dancoff approximation to time-dependent density functional theory within the plane wave basis set framework. *J. Chem. Phys.* **2003**, *118*, 3928–3934.

(41) Curchod, B.; Penfold, T.; Rothlisberger, U.; Tavernelli, I. Nonadiabatic ab initio molecular dynamics using linear-response time-dependent density functional theory. *Open Phys.* **2013**, *11*, 1059–1065.

(42) Mališ, M.; Luber, S. Trajectory Surface Hopping Nonadiabatic Molecular Dynamics with Kohn–Sham ΔSCF for Condensed-Phase Systems. *J. Chem. Theory Comput.* **2020**, *16*, 4071–4086.

(43) Chen, W.-K.; Zhao, X.; Liu, X.-Y.; Xie, X.-Y.; Zeng, Y.; Cui, G. Photoinduced Nonadiabatic Dynamics of a Single-Walled Carbon Nanotube-Porphyrin Complex. *J. Phys. Chem. A* **2024**, *128*, 8803–8815.

(44) Akimov, A. V.; Prezhdo, O. V. The PYXAID Program for Non-Adiabatic Molecular Dynamics in Condensed Matter Systems. *J. Chem. Theory Comput.* **2013**, *9*, 4959–4972.

(45) Smith, B.; Akimov, A. V. A comparative analysis of surface hopping acceptance and decoherence algorithms within the neglect of back-reaction approximation. *J. Chem. Phys.* **2019**, *151*, No. 124107.

(46) Uratani, H.; Nakai, H. Simulating the Coupled Structural–Electronic Dynamics of Photoexcited Lead Iodide Perovskites. *J. Phys. Chem. Lett.* **2020**, *11*, 4448–4455.

(47) Shakiba, M.; Stippell, E.; Li, W.; Akimov, A. V. Nonadiabatic Molecular Dynamics with Extended Density Functional Tight-Binding: Application to Nanocrystals and Periodic Solids. *J. Chem. Theory Comput.* **2022**, *18*, 5157–5180.

(48) Barbatti, M.; Granucci, G.; Persico, M.; Ruckenbauer, M.; Vazdar, M.; Eckert-Maksić, M.; Lischka, H. The on-the-fly surface-hopping program system Newton-X: Application to ab initio simulation of the nonadiabatic photodynamics of benchmark systems. *J. Photochem. Photobiol., A* **2007**, *190*, 228–240.

(49) Barbatti, M.; Ruckenbauer, M.; Plasser, F.; Pittner, J.; Granucci, G.; Persico, M.; Lischka, H. Newton-X: a surface-hopping program for nonadiabatic molecular dynamics. *WIREs Comput. Mol. Sci.* **2014**, *4*, 26–33.

(50) Barbatti, M.; Bondanza, M.; Crespo-Otero, R.; Demoulin, B.; Dral, P. O.; Granucci, G.; Kososki, F.; Lischka, H.; Mennucci, B.; Mukherjee, S.; Pederzoli, M.; Persico, M.; Pinheiro, M.; Jr; Pittner, J.; Plasser, F.; Sangiogo Gil, E.; Stojanovic, L. Newton-X Platform: New Software Developments for Surface Hopping and Nuclear Ensembles. *J. Chem. Theory Comput.* **2022**, *18*, 6851–6865.

(51) Kühne, T. D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V. V.; Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R. Z.; Schütt, O.; Schiffmann, F.; Golze, D.; Wilhelm, J.; Chulkov, S.; Bani-Hashemian, M. H.; Weber, V.; Brštnik, U.; Taillefumier, M.; Jakobovits, A. S.; Lazzaro, A.; Pabst, H.; Müller, T.; Schade, R.; Guidon, M.; Andermatt, S.; Holmberg, N.; Schenter, G. K.; Hehn, A.; Bussy, A.; Belleflamme, F.; Tabacchi, G.; Glöß, A.; Lass, M.; Bethune, I.; Mundy, C. J.; Plessl, C.; Watkins, M.; VandeVondele, J.; Krack, M.; Hutter, J. CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations. *J. Chem. Phys.* **2020**, *152*, No. 194103.

(52) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. *Mol. Phys.* **1997**, *92*, 477–488.

(53) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. *Comput. Phys. Commun.* **2005**, *167*, 103–128.

(54) Iannuzzi, M.; Chassaing, T.; Wallman, T.; Hutter, J. Ground and Excited State Density Functional Calculations with the Gaussian and Augmented-Plane-Wave Method. *CHIMIA* **2005**, *59*, 499.

(55) Strand, J.; Chulkov, S. K.; Watkins, M. B.; Shluger, A. L. First principles calculations of optical properties for oxygen vacancies in binary metal oxides. *J. Chem. Phys.* **2019**, *150*, No. 044702.

(56) Poli, E.; Elliott, J. D.; Chulkov, S. K.; Watkins, M. B.; Teobaldi, G. The Role of Cation-Vacancies for the Electronic and Optical Properties of Aluminosilicate Imogolite Nanotubes: A Non-local, Linear-Response TDDFT Study. *Front. Chem.* **2019**, *7*, No. 210.

(57) Hehn, A.-S.; Sertcan, B.; Belleflamme, F.; Chulkov, S. K.; Watkins, M. B.; Hutter, J. Excited-State Properties for Extended
Systems: Efficient Hybrid Density Functional Methods. *J. Chem. Theory Comput.* **2022**, *18*, 4186-4202.

(58) Sertcan Gökmen, B.; Hutter, J.; Hehn, A.-S. Excited-State Forces with the Gaussian and Augmented Plane Wave Method for the Tamm-Dancoff Approximation of Time-Dependent Density Functional Theory. *J. Chem. Theory Comput.* **2024**, *20*, 8494-8504.

(59) Balasubramani, S. G.; Chen, G. P.; Coriani, S.; Diedenhofen, M.; Frank, M. S.; Franzke, Y. J.; Furche, F.; Grotjahn, R.; Harding, M. E.; Hättig, C.; Hellweg, A.; Helmich-Paris, B.; Holzer, C.; Huniar, U.; Kaupp, M.; Marefat Khah, A.; Karbalaei Khani, S.; Müller, h.; Mack, F.; Nguyen, B. D.; Parker, S. M.; Perlt, E.; Rappoport, D.; Reiter, K.; Roy, S.; Rückert, M.; Schmitz, G.; Sierka, M.; Tapavicza, E.; Tew, D. P.; van Wüllen, C.; Voora, V. K.; Weigend, F.; Wodyński, A.; Yu, J. M. TURBOMOLE: Modular program suite for ab initio quantum-chemical and condensed-matter simulations. *J. Chem. Phys.* **2020**, *152*, 184107.

(60) Franzke, Y. J.; Holzer, C.; Andersen, J. H.; Begušić, T.; Bruder, F.; Coriani, S.; Della Sala, F.; Fabiano, E.; Fedotov, D. A.; Fürst, S.; Gillhuber, S.; Grotjahn, R.; Kaupp, M.; Kehry, M.; Krstić, M.; Mack, F.; Majumdar, S.; Nguyen, B. D.; Parker, S. M.; Pauly, F.; Pausch, A.; Perlt, E.; Phun, G. S.; Rajabi, A.; Rappoport, D.; Samal, B.; Schrader, T.; Sharma, M.; Tapavicza, E.; Treß, R. S.; Voora, V.; Wodyński, A.; Yu, J. M.; Zerulla, B.; Furche, F.; Hättig, C.; Sierka, M.; Tew, D. P.; Weigend, F. TURBOMOLE: Today and Tomorrow. *J. Chem. Theory Comput.* **2023**, *19*, 6859-6890.

(61) Casida, M. E. Time-dependent density-functional theory for molecules and molecular solids. *J. Mol. Struct.:THEOCHEM* **2009**, *914*, 3-18.

(62) Guidon, M.; Hutter, J.; VandeVondele, J. Auxiliary Density Matrix Methods for Hartree-Fock Exchange Calculations. *J. Chem. Theory Comput.* **2010**, *6*, 2348-2364.

(63) Mitrić, R.; Werner, U.; Bonačić-Koutecký, V. Nonadiabatic dynamics and simulation of time resolved photoelectron spectra within time-dependent density functional theory: Ultrafast photo-switching in benzylideneaniline. *J. Chem. Phys.* **2008**, *129*, 164118.

(64) Werner, U.; Mitrić, R.; Suzuki, T.; Bonačić-Koutecký, V. Nonadiabatic dynamics within the time dependent density functional theory: Ultrafast photodynamics in pyrazine. *Chem. Phys.* **2008**, *349*, 319-324. Electron Correlation and Molecular Dynamics for Excited States and Photochemistry.

(65) Mitrić, R.; Werner, U.; Wohlgemuth, M.; Seifert, G.; Bonačić-Koutecký, V. Nonadiabatic Dynamics within Time-Dependent Density Functional Tight Binding Method. *J. Phys. Chem. A* **2009**, *113*, 12700-12705.

(66) Plasser, F.; Crespo-Otero, R.; Pederzoli, M.; Pittner, J.; Lischka, H.; Barbatti, M. Surface Hopping Dynamics with Correlated Single-Reference Methods: 9H-Adenine as a Case Study. *J. Chem. Theory Comput.* **2014**, *10*, 1395-1405.

(67) Grimme, S.; Bannwarth, C.; Shushkov, P. A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1 - 86). *J. Chem. Theory Comput.* **2017**, *13*, 1989-2009.

(68) Bannwarth, C.; Caldeweyher, E.; Ehlert, S.; Hansen, A.; Pracht, P.; Seibert, J.; Spicher, S.; Grimme, S. Extended tight-binding quantum chemistry methods. *WIREs Comput. Mol. Sci.* **2021**, *11*, No. e01493.

(69) Iannuzzi, M.; Wilhelm, J.; Stein, F.; Bussy, A.; Elgabarty, H.; Golze, D.; Hehn, A.; Graml, M.; Marek, S.; Gökmen, B. S.; Schran, C.; Forbert, H.; Khaliullin, R. Z.; Kozhevnikov, A.; Taillefumier, M.; Meli, R.; Rybkin, V.; Brehm, M.; Schade, R.; Schütt, O.; Pototschnig, J. V.; Mirhosseini, H.; Knüpfer, A.; Marx, D.; Krack, M.; Hutter, J.; Kühne, T. D. CP2K Program Package Made Simple 2025 https://arxiv.org/abs/2508.15559.

(70) Tomasello, G.; Humeniuk, A.; Mitrić, R. Exploring Ultrafast Dynamics of Pyrazine by Time-Resolved Photoelectron Imaging. *J. Phys. Chem. A* **2014**, *118*, 8437-8445.

(71) Stock, G.; Domcke, W. Detection of ultrafast molecular-excited-state dynamics with time- and frequency-resolved pump-probe spectroscopy. *Phys. Rev. A* **1992**, *45*, 3032-3040.

(72) Richings, G. W.; Habershon, S. MCTDH on-the-fly: Efficient grid-based quantum dynamics without pre-computed potential energy surfaces. *J. Chem. Phys.* **2018**, *148*, No. 134116.

(73) Belyaev, A. K.; Lasser, C.; Trigila, G. Landau-Zener type surface hopping algorithms. *J. Chem. Phys.* **2014**, *140*, No. 224108.

(74) Xie, W.; Sapunar, M.; Došlić, N.; Sala, M.; Domcke, W. Assessing the performance of trajectory surface hopping methods: Ultrafast internal conversion in pyrazine. *J. Chem. Phys.* **2019**, *150*, No. 154119.

(75) Krack, M.; Parrinello, M. All-electron ab-initio molecular dynamics. *Phys. Chem. Chem. Phys.* **2000**, *2*, 2105-2112.

(76) Becke, A. D. A new mixing of Hartree-Fock and local density-functional theories. *J. Chem. Phys.* **1993**, *98*, 1372-1377.

(77) Weigend, F.; Furche, F.; Ahlrichs, R. Gaussian basis sets of quadruple zeta valence quality for atoms H-Kr. *J. Chem. Phys.* **2003**, *119*, 12753-12762.

(78) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. *Phys. Chem. Chem. Phys.* **2005**, *7*, 3297.

(79) Lippert, G.; Hutter, J.; Parrinello, M. The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations. *Theor. Chem. Acc.* **1999**, *103*, 124-140.

(80) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. *J. Chem. Phys.* **2007**, *127*, No. 114105.

(81) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. *Phys. Rev. B* **1996**, *54*, 1703-1710.

(82) Crespo-Otero, R.; Barbatti, M. Spectrum simulation and decomposition with nuclear ensemble: formal derivation and application to benzene, furan and 2-phenylfuran. *Theor. Chem. Acc.* **2012**, *131*, 1237.

(83) Github. Computational data. https://github.com/SPECATgroup/NEWTONXinterface.

(84) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. *Phys. Rev. Lett.* **1996**, *77*, 3865-3868.

(85) Brémond, E.; Savarese, M.; Adamo, C.; Jacquemin, D. Accuracy of TD-DFT Geometries: A Fresh Look. *J. Chem. Theory Comput.* **2018**, *14*, 3715-3727.

(86) Kretz, B.; Egger, D. A. Accurate non-adiabatic couplings from optimally tuned range-separated hybrid functionals. *J. Chem. Phys.* **2022**, *157*, No. 101104.

(87) Suzuki, Y.-I.; Fuji, T.; Horio, T.; Suzuki, T. Time-resolved photoelectron imaging of ultrafast S2→S1 internal conversion through conical intersection in pyrazine. *J. Chem. Phys.* **2010**, *132*, No. 174302.

(88) Stert, V.; Farmanara, P.; Radloff, W. Electron configuration changes in excited pyrazine molecules analyzed by femtosecond time-resolved photoelectron spectroscopy. *J. Chem. Phys.* **2000**, *112*, 4460-4464.

(89) Kanno, M.; Ito, Y.; Shimakura, N.; Koseki, S.; Kono, H.; Fujimura, Y. Ab initio quantum dynamical analysis of ultrafast nonradiative transitions via conical intersections in pyrazine. *Phys. Chem. Chem. Phys.* **2015**, *17*, 2012-2024.

(90) Yamazaki, I.; Murao, T.; Yamanaka, T.; Yoshihara, K. Intramolecular electronic relaxation and photoisomerization processes in the isolated azabenzene molecules pyridine, pyrazine and pyrimidine. *Faraday Discuss. Chem. Soc.* **1983**, *75*, 395-405.

(91) Azumi, T. Defect Phosphorescence of Pyrazine Crystal. II. The Nature of the Defect. *J. Chem. Phys.* **1971**, *54*, 85-91.

(92) Moomaw, W. R.; El-Sayed, M. A. Phosphorescence of Crystalline Pyrazine at 4.2°K. *J. Chem. Phys.* **1968**, *48*, 2502-2514.

(93) Angulo, G.; Dobkowski, J.; Kapturkiewicz, A. Luminescence properties of diamino-dicyano substituted benzene and 1,4-pyrazine. *J. Photochem. Photobiol., A* **2011**, *225*, S2-S7.
(94) Shirai, K.; Yanagisawa, A.; Takahashi, H.; Fukunishi, K.; Matsuoka, M. Syntheses and fluorescent properties of 2,5-diamino-3,6-dicyanopyrazine dyes. *Dyes Pigm.* **1998**, *39*, 49–68.

(95) Hirosawa, R.; Akune, Y.; Endo, N.; Hatano, S.; Hosokai, T.; Sato, H.; Matsumoto, S. A variety of solid-state fluorescence properties of pyrazine dyes depending on terminal substituents. *Dyes Pigm.* **2017**, *146*, 576–581.

(96) Rajamalli, P.; Chen, D.; Suresh, S. M.; Tsuchiya, Y.; Adachi, C.; Zysman-Colman, E. Planar and Rigid Pyrazine-Based TADF Emitter for Deep Blue Bright Organic Light-Emitting Diodes. *Eur. J. Org. Chem.* **2021**, *2285–2293*.

(97) Singh, P. S.; Chacko, S.; Kamble, R. M. Dicyanopyrazino phenanthrene based charge transfer derivatives: Role of amine donor in tuning of photophysical, aggregation-induced emission, electrochemical and theoretical properties. *J. Mol. Struct.* **2023**, *1271*, No. 134052.

(98) Lewis, T. P.; Ragin, H. R. Polarized single-crystal absorption spectra of pyrazine and tetramethylpyrazine. *J. Am. Chem. Soc.* **1972**, *94*, 5566–5573.

(99) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. *Phys. Rev. B* **1998**, *58*, 3641–3662.

(100) Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules. *J. Chem. Phys.* **2012**, *137*, No. 014512.

<figure><img src="image_16.png" /></figure>