# COLUMBUS An Efficient and General Program Package for Ground and Excited State Computations Including Spin-Orbit Couplings and Dynamics

Felix Plasser, Hans Lischka, Ron Shepard, Péter G Szalay, Russell M Pitzer, Rodolpho L R Alves, Adelia J A Aquino, Jochen Autschbach, Mario Barbatti, Jhonatas R Carvalho, et al.

► To cite this version:

Felix Plasser, Hans Lischka, Ron Shepard, Péter G Szalay, Russell M Pitzer, et al. COLUMBUS An Efficient and General Program Package for Ground and Excited State Computations Including Spin-Orbit Couplings and Dynamics. Journal of Physical Chemistry A, 2025, 10.1021/acs.jpca.5c02047. [hal-05153802](https://hal.science/hal-05153802v1)

HAL Id: hal-05153802

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

Submitted on 9 Jul 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>

# COLUMBUS—An Efficient and General Program Package for Ground and Excited State Computations Including Spin-Orbit Couplings and Dynamics

Published as part of *The Journal of Physical Chemistry A* **special issue** "*Quantum Chemistry Software for Molecules and Materials*".

Felix Plasser,* Hans Lischka,* Ron Shepard, Péter G. Szalay, Russell M. Pitzer, Rodolpho L. R. Alves, Adelia J. A. Aquino, Jochen Autschbach, Mario Barbatti, Jhonatas R. Carvalho, Julio C. V. Chagas, Leticia González, Andreas Hansen, Bhumika Jayee, Miklos Kertesz, Francisco B. C. Machado, Spiridoula Matsika, Silmar A. do Monte, Saikat Mukherjee, Dana Nachtigallová, Reed Nieman, Vytor P. Oliveira, Markus Oppel, Carol A. Parish, Jiri Pittner, Luan G. F. dos Santos, Armin Scrinzi, Mahesh K. Sit, Rene F. K. Spada, Mushir Thodika, Daniel C. A. Valente, Álvaro Vázquez-Mayagoitia, Elizete Ventura, Julia Westermayr, Aleksandr Zaichenko, and Zhiyong Zhang

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

Cite This: https://doi.org/10.1021/acs.jpca.5c02047

<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>

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

Article Recommendations

**ABSTRACT:** The COLUMBUS program system provides the tools for performing high-level multireference (MR) computations, including the multireference configuration interaction (MRCI) method and its multireference averaged quadratic coupled cluster (MR-AQCC) extension, allowing computations on a wide range of fascinating atomic and molecular systems, including the treatment of open-shells and complicated excited state phenomena. The inclusion of spin-orbit coupling (SOC) directly within the MRCI step enables the description of systems containing heavy elements, such as lanthanides and actinides, whose properties are strongly influenced by SOC. Analytic energy gradients and nonadiabatic couplings at the correlated MRCI level provide the foundation for a variety of dynamics studies, giving insight into ultrafast photochemistry. New and ongoing method developments in COLUMBUS include the computation of spin densities, improved descriptions of ionic states, enhancements to the AQCC method, and the porting of COLUMBUS to graphical processing units (GPUs). New external interfaces enable an enhanced description of electronic resonances and molecules in strong laser fields. This work highlights these new developments while providing a detailed account of the diverse applications of COLUMBUS in recent years.

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

## 1. INTRODUCTION

Within the last decades, quantum chemistry has become an indispensable component of molecular science investigations. Whereas ground-state reactivity can often be described using standard density functional theory (DFT) computations, there are significant and vital classes of problems where DFT fails, most notably in cases involving strongly correlated electrons, distorted structures, broken bonds, and doubly excited states.<sup>1,2</sup> In such cases, multireference (MR) methods provide an attractive option by allowing the treatment of all relevant electronic configurations of a given molecular system on an equal footing. Aside from accurate energies, it is often also essential to be able to compute other electronic properties,

including energy gradients, allowing one to determine molecular geometries and their dynamic evolution, transition dipole moments for optical transition strengths, as well as nonadiabatic (NAC) and spin-orbit coupling (SOC) matrix elements to determine interstate transition probabilities.

Received: March 27, 2025

Revised: June 23, 2025

Accepted: June 24, 2025

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 1.** Summary of the functionality provided within COLUMBUS. Features in the release version are marked with solid lines, while features under development are marked with dashed lines.</figcaption></figure>

In the spirit of the above discussion, the purpose of COLUMBUS is to provide a platform for flexible and high-level multireference computations along with the functionality to compute the relevant properties and interstate coupling terms. Within the landscape of multireference programs, COLUMBUS stands out through its flexible multiconfigurational self-consistent field (MCSCF) module and its flexible, powerful, and well-parallelized multireference configuration interaction (MRCI) module. The latter includes an extension for variational spin-orbit configuration interaction and various extensivity corrections. Furthermore, COLUMBUS provides energy gradients, nonadiabatic couplings, spin-orbit couplings, and other electronic properties at the MRCI and MCSCF levels.

The purpose of this paper is to present the computational capabilities of COLUMBUS, including ongoing developments, and to illustrate them through a range of recent applications. We start with a brief overview of the program structure and main capabilities of Columbus (Section 2). Subsequently, Section 3 introduces recent advances, including the computation of MCSCF spin densities within the originally spin-free graphical unitary group approach (GUGA), new treatments of ionic states as defined in valence-bond theory—long recognized as a challenge³⁻⁵ for many multireference frameworks—and a detailed account of the averaged quadratic coupled cluster (AQCC) approach within a linear response formalism (AQCC-LRT). This section also covers current efforts to port COLUMBUS to GPUs and highlights two new external interfaces, used for computing electronic resonances and modeling attosecond experiments, respectively. Section 4 presents selected case studies where multireference methods are essential, including calibration of the semiempirical fractional occupation number weighted electron density (FOD)⁶ method against multireference unpaired densities and applications to graphene defects, phenalenyl oligomers, polyradicals, and ionic hydrocarbons. We conclude this part

with the intricate case of a uranium-uranium bond, where both multireference and SOC effects are significant. Section 5 focuses on excited-state treatments, beginning with the analysis of covalent and ionic character and continuing with studies on doped periacenes and the triatomic CUO molecule, again emphasizing spin-orbit coupling effects. Finally, Section 6 demonstrates the use of COLUMBUS in various photo-dynamics simulations, including recent developments that leverage trajectory data for machine learning applications.

## 2. OVERVIEW

A summary of the COLUMBUS program structure is given in Figure 1. We describe these features here and refer to the various sections that exemplify their use. At the heart of COLUMBUS are its MCSCF and MRCI modules. Both are based on the GUGA formalism,⁷ which enables a flexible definition of various wave function types and an efficient computation of the associated matrix elements. These methods enable the accurate treatment of challenging electronic structures, such as those found in graphene defects (Section 4.2), polyradicals (Section 4.3), and nonstandard bonding situations (Section 4.5). In MCSCF and MRCI computations, COLUMBUS allows for complete active space SCF (CASSCF) and restricted active space SCF (RASSCF), as well as more specialized and highly customized active space definitions (one such example is described in Section 6.1). COLUMBUS utilizes an uncontracted MRCI with singles and doubles (MR-CISD) approach⁸,⁹ implemented using a direct configuration interaction (CI) algorithm, which avoids the explicit calculation of the Hamiltonian matrix elements. Combined with an efficient parallelization algorithm, CI dimensions of several billion can be handled routinely on standard computer clusters, making many interesting calculations possible,¹⁰⁻¹² as illustrated in various examples below. Further improvement, aimed at utilizing GPUs, is currently ongoing (Section 3.4).
A significant advantage of uncontracted MRCI is that, due to its variational nature, it is relatively straightforward to compute electronic (transition) properties, gradients, and nonadiabatic couplings (shown to the right in Figure 1), all of which are implemented for MCSCF and MRCI. The availability of analytic gradients allows for structure optimizations (illustrated in various sections below), and the combination with nonadiabatic couplings provides the basis for photodynamics simulations (Section 6). Along with the wave function properties, COLUMBUS naturally produces natural orbitals, densities, transition densities, and so on for the correlated wave functions. Recently, we have also demonstrated how the originally spin-free GUGA formalism can nonetheless be utilized to compute spin-densities (Section 3.1).

COLUMBUS features functionality to treat relativistic spin-orbit coupling (SOC) effects within its SO-CI module, also based on the GUGA formalism.¹³ SO-CI can be performed in a variational fashion, meaning that SOC, static, and dynamic correlation are treated in one combined step, providing a high-level description even of challenging electronic structures. Using SO-CI becomes particularly important for heavy elements with unpaired electrons, and its applications are illustrated below in the case of bonding between two uranium atoms (Section 4.6) and the CUO triatomic molecule (Section 5.3). The computation of gradients and nonadiabatic couplings from SO-CI is available in a development version but has not been officially released at the time of writing. Alternatively, a perturbative version of the SO-CI module is available,¹⁴ where SOC is computed between standard spin-free MRCI wave functions. This implementation is particularly well-suited to dynamics simulations, as presented in Section 6.2.

CI methods are advantageous due to their variational nature, and this is reflected in numerically stable results as well as the relatively straightforward implementation of gradients, non-adiabatic couplings, and other wave function properties. Conversely, a difficulty when applying CI is its lack of extensivity, meaning that electron correlation in larger molecules is not as well described as in smaller molecules. Practically, one finds that MRCI is generally reliable only if an extensivity correction is used and COLUMBUS provides various possible options (Figure 1, upper left). As one conceptually simple and computationally cheap option, a *posteriori* corrections can be used.¹⁵,¹⁶ These are often generically denoted as (+Q). Within COLUMBUS, different variations are implemented, and we generally recommend¹⁷ using the Pople (+P)¹⁸ or Davidson-Silver (+DV3)¹⁹ corrections. A *posteriori* corrections tend to provide reliable energies but it is not possible to use them to compute energy gradients. For this reason, COLUMBUS also provides iterative extensivity corrections, that is, methods where already the CI procedure is modified. The most prominent one of these is the MR-averaged quadratic coupled cluster (MR-AQCC) method.²⁰ MR-AQCC provides reliable energies of the ground and excited states along with their gradients. The MR-AQCC method and its linear response theory (LRT) variant are discussed in Section 3.3. Furthermore, the practical influence of extensivity corrections in producing accurate excitation energies is discussed in detail in Sections 3.2 and 5.1 within a valence-bond framework focusing on the ionic and covalent character of the state. Finally, noting that MR-AQCC gives access to extensivity-corrected density matrices, we use it in Section 4.1 to provide an accurate reference for calibrating an approximate method for describing unpaired electrons.

One crucial structural feature of COLUMBUS is the modular concept, which organizes different computational steps, including atomic orbital (AO) integrals, GUGA and wave function setup, integral transformation, CI steps, and analytic energy gradient and nonadiabatic coupling calculations. This approach enables the combination of COLUMBUS with other quantum chemical program packages to share atomic orbital (AO) integrals. In addition to the proprietary Argos program, well-established interfaces to DALTON²¹ and OpenMolcas²² are available, along with a more experimental interface to NWChem²³ (see Figure 1, bottom). Furthermore, the SCF, MCSCF, and geometry optimization modules of OpenMolcas can be used in combination with the MRCI program module of COLUMBUS, providing enhanced flexibility for setting up various jobs.²⁴ The modular setup also provides an ideal basis for external interfaces. Below, we illustrate this in the case of complex absorbing potentials used for electronic resonances (Section 3.5) and hybrid antisymmetrized coupled channels for attosecond spectroscopy (Section 3.6).

COLUMBUS is developed as an open-source project via a public repository²⁵ and is released under the GNU Lesser General Public License. Releases are distributed from the same web page²⁵ in the form of source code and precompiled executables. COLUMBUS is primarily developed for Linux systems. It has been installed and used on a variety of high-performance cluster architectures.

### 3. METHOD DEVELOPMENTS

#### 3.1. Spin Densities at the MCSCF Level.
Recently, the calculation of spin densities was implemented in COLUMBUS for wave functions obtained via MCSCF calculations.²⁶ Given a wave function $|\psi; S, M\rangle$, which is an eigenfunction of the operator $\hat{S}^2$, that is

$$ \hat{S}^2 |\psi; S, M\rangle = S(S+1)|\psi; S, M\rangle \quad (1) $$

with eigenvalues $S \ge 0$ and of the operator $\hat{S}_z$

$$ \hat{S}_z |\psi; S, M\rangle = M|\psi; S, M\rangle \quad (2) $$

determining the projection along the z-axis with quantum numbers $M = -S, -S+1, \dots, S-1, S$, the calculation of spin-density matrix elements is given by

$$ D_{qp}^{(1,0,M)} = \langle \psi; S, M | a_{p\alpha}^{\dagger} a_{q\alpha} - a_{p\beta}^{\dagger} a_{q\beta} | \psi; S, M \rangle \quad (3) $$

In COLUMBUS, this calculation is implemented using the Graphical Unitary Group Approach (GUGA),²⁷ as this procedure avoids transforming the wave function from a configuration state function (CSFs) basis to a determinantal basis to carry out the difference of the $\alpha$ and $\beta$ blocks. Note that this last procedure would be highly inefficient or unfeasible for large CSF expansions since the determinantal expansion grows faster than the CSF one, especially for larger numbers of open-shell electrons.

The current implementation calculates the spin-density matrix elements for $M = S$ using²⁶,²⁸

$$ D_{qp}^{(1,0,S)} = \frac{(2 - \frac{1}{2}N)}{(S + 1)} D_{qp}^{(0)} - \frac{1}{(S + 1)} \sum_{k} d_{qkkp} \quad (4) $$

in which $N$ is the total number of electrons, $D_{qp}^{(0)}$ is an element from the one-particle reduced density matrix (1-RDM) and $d_{qkkp} = \langle \psi | \hat{e}_{pkkq} | \psi \rangle$ are the unsymmetrized two-particle reduced
density matrix (2-RDM) elements. Both the 1-RDM and 2-RDM are available within the GUGA formalism. Note that the calculation of the last term on the right-hand side of eq 4 would require a sum over elements of a four-dimensional array. However, a more efficient approach was implemented to compute the spin-density matrix alongside the 1-RDM and 2-RDM to reduce the additional computational effort to obtain the elements $D_{qp}^{(1,0,S)}$.

The analysis of the spin density is useful, for example, for the study of organic radicals. In particular, diradicals have applications in various areas of optoelectronics and spintronics.²⁹⁻³³ In this work, we employed the implementation to calculate the spin density for the diradical para-quinodimethane. Investigations into the electronic structure of this system can be traced to the works by Szwarc³⁴ and Coulson et al.³⁵ Recent work by Chagas et al. has also explored the electronic effects of the CH₂ groups on the benzenoid ring.³⁶ A detailed study of the electronic structure of this molecule was recently carried out by Matasović et al.³⁷ to investigate the conversion from a closed-shell character attained with the CH₂ groups coplanar to the ring structure to an open-shell character with the CH₂ groups orthogonal. Here, we utilized the planar geometry from that work, optimized at the PBE/ANO-S-VDZP level, and two additional structures with different values for the θ angle between the ring plane and the CH₂ group (Figure 2a): 0, 50, and 90°.

<figure><img src="image_11.png" /><figcaption>**Figure 2.** para-Quinodimethane molecule with panel (a) highlighting the angle (θ) between the ring plane and the CH₂ group. Panels (b–d) display the spin densities for the T₁ state of the studied structures with different θ angles. The employed isovalue is equal to 0.003 eÅ⁻³.</figcaption></figure>

Next, we carried out calculations using the SA-CASSCF method, considering 8 electrons in 8 orbitals (CAS(8,8)) and the cc-pVDZ basis set for the lowest triplet state (T₁). The spin densities for the three structures are presented in Figure 2.

From the inspection of Figure 2, the decrease in the spin density of the ring structure is noticeable with the increase in the angle θ. This highlights the loss in delocalization upon twisting and illustrates that the spin density is increasingly localized on the radical centers.

## 3.2. Diagnostics and Scaling Procedures for Ionic States.

Whereas MCSCF presents itself as a compelling method for providing qualitatively correct starting points in various situations, including open shells and highly excited states, one critical area that poses notorious challenges is the description of ionic $ππ*$ states.³⁻⁵ These states can be best understood in a valence-bond framework.³⁸,³⁹ To summarize, within the standard MO picture, the singlet and triplet highest-occupied molecular orbital/lowest-unoccupied molecular orbital (HOMO/LUMO) excited states are given as

$$ \Psi_s = |h\bar{l}\rangle + |l\bar{h}\rangle \quad (5) $$

$$ \Psi_T = |h\bar{l}\rangle - |l\bar{h}\rangle \quad (6) $$

where h and l refer to the HOMO and LUMO, and the bar signifies β-spin. Using a valence-bond-type description and expressing these equations in terms of localized orbitals, denoted a and b, yields¹⁷,³⁹

$$ \Psi_s = |a\bar{a}\rangle - |b\bar{b}\rangle \quad (7) $$

$$ \Psi_T = |b\bar{a}\rangle - |a\bar{b}\rangle \quad (8) $$

Crucially, the singlet consists of configurations in which both electrons occupy either orbital a or orbital b, and the state is therefore classified as ionic. The opposite is true for the triplet, which is classified as a diradical. Practically speaking, the ionic nature of the singlet requires enhanced electron correlation in its description, often posing severe challenges to the MCSCF method. There are two interrelated problems: how do we diagnose issues related to ionic states, and how do we fix these problems? Both have been addressed recently, as discussed below.

A new diagnostic for ionic states, denoted $Q_a^t$, has been developed in ref 17. $Q_a^t$ is based on a population analysis of the transition density, carried out with TheoDORE,⁴⁰ to postprocess the COLUMBUS computation. Before proceeding, we comment on the range of values accessible by $Q_a^t$. It is per construction always greater than zero. However, there is no obvious upper limit. Whereas $Q_a^t$ values of singlet states are usually found to be below one, one finds values close to two for some triplet states.¹⁷

**Figure 3** presents an analysis of computed errors in excitation energies (using the theoretical best estimates from ref 41 as reference) for a test set of 11 π-conjugated molecules plotted against the $Q_a^t$ diagnostic. Starting with the SA-CASSCF level (Figure 3a) and considering the $ππ*$ states (red circles), we find that the diagnostic performs precisely as intended. For the covalent $ππ*$ states seen at low $Q_a^t$ values (below 0.2), we find good accuracy, and all errors are below 0.5 eV. Conversely, for the ionic $ππ*$ states seen at higher $Q_a^t$ values (above 0.3), significant errors of 1–2 eV are observed. Finally, the $nπ*$ states (green squares) yield intermediate errors. Figure 3a thus highlights the power of the $Q_a^t$ diagnostic for $ππ*$ states, providing a clear differentiation between cases
<figure><img src="image_12.png" /><figcaption>**Figure 3.** Errors of computed vertical singlet excitation energies for a test set of 11 molecules plotted against the $Q_a^t$ diagnostic measuring ionic character using (a) SA-CASSCF, (b) MR-CIS, (c) MR-CISD, and (d) MR-CISD+P, all using the aug-cc-pVDZ basis set. States are grouped according to type: $\pi\pi^*$ (circles) and $n\pi^*$ (squares) states; the out-of-plane $\pi\pi^*$ state of cyanoformaldehyde is shown as an empty circle. Adapted from ref 17. Copyright 2023 American Chemical Society.</figcaption></figure>

where the SA-CASSCF is expected to be reliable and where it is not.

Progressing along the hierarchy of MRCI methods, as implemented in COLUMBUS, generally improves the results from CASSCF, and this is shown in Figure 3b-d. However, an important point to note is that MR-CIS sometimes yields improved results compared to MR-CISD, and this discrepancy can be attributed to the lack of size extensivity. But, using finally an extensivity-corrected version (MR-CISD+P) yields excellent results for all types of states.

While the above work demonstrates that ionic states can be effectively detected, correcting their energies at the CASSCF stage would be beneficial. Work in this direction is currently in progress. The formalism starts by realizing that the energies of the singlet and triplet states in eqs 5 and 6 differ by twice the exchange integral $g_{hl,hl}$ where the triplet is stabilized and the singlet is raised in energy.42,43 Furthermore, it is worth noting that $\sigma$-correlation can be understood to act by screening this exchange interaction.17,43,44 The idea behind the new procedure is to scale the involved terms to mimic screening via the $\sigma$-electrons. Practically, we can define the diagonal shift Hamiltonian

$$ \hat{H}_{\text{shift}} = (1 - \mu) \sum_{p<q} \hat{p}_{pq} g_{pq, pq} \quad (9) $$

where $\hat{p}_{pq}$ is the projector onto all CSFs where orbitals $p$ and $q$ are singly occupied and singlet spin-coupled, and all other orbitals are either unoccupied or doubly occupied. The value of $\mu$ is left as an adjustable parameter. A scaled MCSCF version can now be defined by subtracting this shift Hamiltonian from the usual MCSCF Hamiltonian. Preliminary results suggest that this procedure specifically targets singlet $\pi\pi^*$ states, leaving the ground state and other excited states largely unaffected.

The proposed procedure directly targets the underlying two-electron integrals and is expected to allow for a much more fine-grained tuning than global energy scaling procedures.45,46 Crucially, as opposed to multiconfigurational pair density functional theory,47 the shift is included self-consistently in the MCSCF procedure, meaning that it also affects orbital optimization and ultimately enables the computation of energy gradients and nonadiabatic couplings. Work on this is currently in progress, and the developed method will be presented in a separate communication.

### 3.3. AQCC-LRT as a Viable Alternative to MR-CISD Calculations of Excited States. Multireference Averaged Coupled Pair Functional (MR-ACPF)48 and MR-AQCC20,49

methods are powerful tools for correcting the size-extensivity error of MR-CISD. They have been utilized in numerous applications (for a review, see refs 50,51). Their main advantage over the alternative, Davidson-type corrections, is that a density is associated with the corrected energy,48 and therefore, gradient and other properties can be calculated.52,53

In the MR-AQCC method,20,49 a functional of the correlation energy ($\Delta E_\alpha$) of state $\alpha$

$$ F_{\alpha}(\epsilon^{\alpha}) = \frac{\langle \sum_{i} c_{i}^{\alpha} \Phi_{i} | \hat{H} - E_{0}^{\alpha} | \sum_{i} c_{i}^{\alpha} \Phi_{i} \rangle}{\sum_{i \in \text{int}} (c_{i}^{\alpha})^{2} + G \sum_{i \notin \text{int}} (c_{i}^{\alpha})^{2}} \quad (10) $$

is optimized to get the coefficients ($c_i^\alpha$) corresponding to the expansion functions $\Phi_i$. Here $E_0^\alpha$ is the reference energy of the state considered, and the two summations in the denominator run for the internal and all other excitations, respectively. By choosing the factor $G$ appropriately, various methods can be defined. Trivially, $G=1$ gives the MR-CISD functional. For MR-AQCC, we have $G=1 - \frac{(n_e-3)(n_e-2)}{n_e(n_e-1)}$, while for MR-ACPF $G = \frac{2}{n_e}$ with $n_e$ being the number of (correlated) electrons.

The stationary equations of this functional can be cast into a matrix eigenvalue problem54,55

$$ \langle \Phi_j | (\hat{H} - E_0^\alpha + \hat{\Delta}_\alpha) | \sum_i c_i^\alpha \Phi_i \rangle = \Delta E_\alpha c_j^\alpha \quad (11) $$

with the diagonal shift being

$$ \hat{\Delta}_{\alpha} = \sum_{k \notin \text{int}} (1 - G) \Delta E_{\alpha} |\Phi_k\rangle \langle \Phi_k| \quad (12) $$

Note that for CI, the shift is zero. Thus, the formalism reproduces the usual CI eigenvalue equations. On the other hand, for $G \neq 1$, the shift is state-dependent since it includes the correlation energy of the given state. In this case, the eigenvalues ($\Delta E_\alpha$) and the coefficient vectors ($\epsilon^\alpha$) belonging to different states are eigenpairs of different Hamiltonians. Whereas the resulting nonorthogonality of the eigenvectors is not a problem for the calculation of single-state properties,48 it does not allow for the computation of transition properties.

The MR-AQCC method has been successfully applied to numerous ground-state applications, such as describing the biradical character and singlet-triplet splitting in polycyclic aromatic hydrocarbons (PAHs)11,50,57 as well as the study of bonding in the chromium dimer.10 It can also be applied to excited states58 without symmetry or spin restriction. This is achieved in COLUMBUS by a root-following technique, where
<table><caption>Table 1. Relative Energy (eV), Transition Dipole Moments (μ^t, a.u.) and Q_u^t Diagnostic for the Three Low-Lying Excited States of All-trans Hexatriene Obtained by the Different Methods</caption><thead><tr><th rowspan="2">state</th><th colspan="3">MR-CISD</th><th colspan="2">MR-CISD+P</th><th colspan="2">MR-AQCC</th><th colspan="2">MR-AQCC-LRT</th></tr><tr><th>ΔE</th><th>μ<sup>t</sup></th><th>Q<sub>u</sub><sup>t</sup></th><th>ΔE</th><th>ΔE</th><th>ΔE</th><th>ΔE</th><th>μ<sup>t</sup></th><th>Q<sub>u</sub><sup>t</sup></th></tr></thead><tbody><tr><td>2<sup>1</sup>A<sub>g</sub></td><td>5.816</td><td>0.000</td><td>0.088</td><td>5.755</td><td>5.743</td><td>5.755</td><td>0.000</td><td>0.077</td><td></td></tr><tr><td>1<sup>1</sup>B<sub>u</sub></td><td>6.520</td><td>2.826</td><td>0.828</td><td>5.802</td><td>5.764</td><td>6.197</td><td>2.745</td><td>0.789</td><td></td></tr><tr><td>2<sup>1</sup>B<sub>u</sub></td><td>7.005</td><td>0.129</td><td>0.131</td><td>6.915</td><td>6.880</td><td>6.902</td><td>0.145</td><td>0.071</td><td></td></tr></tbody></table>

the subspace root with the maximal overlap with a selected reference root is updated during the Davidson subspace approach.⁵⁵ This procedure works well when the desired roots are well-separated and possess distinct characteristics. However, in practice, it lacks the stability of MR-CISD calculations, and intruders often hinder convergence.

We adopted a modified version of the AQCC method, which had been previously developed,⁵⁴ and removes the disadvantage of AQCC associated with the state-dependent diagonal shift. However, it still has the advantage of including size-extensivity corrections intrinsically. As shown in ref 54, the state dependence of the diagonal shift can be avoided by an alternative derivation of the excited state equations using linear response theory (LRT).⁵⁹ Starting from the AQCC ground state solution, for the excited states, the eigenvalue equations of the response matrix⁵⁹ should be considered, i.e., the following equations need to be solved

$$ \langle \Phi_i | (\hat{H} - E_0^\alpha + \hat{\Delta}_0) | \sum_i c_i^\alpha \Phi_i \rangle = \Delta E_\alpha c_j^\alpha \quad (13) $$

for more details on the derivation of the method called AQCC-LRT, see ref 54 here, we draw attention only to the fact that for all excited states, the shift depends on the correlation energy of the ground state (ΔE₀), i.e., energies and coefficients of all excited states are obtained as the eigenpairs of the same matrix. One can, therefore, calculate excited states just like from CI and obtain several excited states at the same time without the necessity of using root following.

AQCC and AQCC-LRT would give the same results if the reference energies for the ground and excited states were of the same quality, resulting in (nearly) identical correlation energy. Although this is often not the case, ref 54 shows that the excitation energies are usually very similar to the two AQCC procedures.

The other advantage of AQCC-LRT is that using the LRT technology,⁵⁹ it is possible to define the transition density between state 0 and α as⁵⁴

$$ \Gamma_t^{0\alpha} = a_0^{\text{CI}} \Gamma_t^{0\alpha} + (1 - a_0)^{\text{ref}} \Gamma_t^{0\alpha} \quad (14) $$

with

$$ a_0 = \frac{\sum_i (c_i^0)^2}{\sum_{i \in \text{int}} (c_i^0)^2 + G \sum_{i \notin \text{int}} (c_i^0)^2} \quad (15) $$

where $^{\text{CI}}\Gamma_t^{0\alpha}$ and $^{\text{ref}}\Gamma_t^{0\alpha}$ are the CISD and reference transition densities, respectively.

Compared with standard MR-AQCC calculations, below, we demonstrate the usefulness of MR-AQCC-LRT by studying all-trans hexatriene where the different character of the lowest excited states poses a particular challenge. The condition for the MR-AQCC and MR-AQCC-LRT to yield similar results is not generally fulfilled, as the description of the ionic state in the reference space is significantly worse; thus, the correlation

energy is expected to be larger in magnitude than for the ground state. In addition, we demonstrate that the transition properties can also be calculated with the MR-AQCC-LRT method, allowing a clear characterization of the wave function.

The ground and the three lowest excited states (2<sup>1</sup>A<sub>g</sub>, 1<sup>1</sup>B<sub>u</sub> and 2<sup>1</sup>B<sub>u</sub>) of the all-trans hexatriene are considered below. Reported results are from MR-CISD, MR-CISD with Pople's approximate extensivity correction (MR-CISD+P),¹⁸ MR-AQCC, and MR-AQCC-LRT. The geometries and orbitals are taken from the calculations reported below (see Section 5.1). The reference space was always a full π-valence CAS(6,6), and only the core orbitals were frozen. All the calculations used the cc-pVDZ basis set.⁶⁰ For the analysis of the transition properties, the TheODORE 3.2 program⁴⁰ was used. We present the transition densities and the Q<sub>u</sub><sup>t</sup> diagnostic,¹⁷ as also discussed in the previous section. In this context it is worth noting that Q<sub>u</sub><sup>t</sup> values and transition dipole moments can only be computed using MR-CISD and MR-AQCC-LRT while this is not possible with MR-CISD+P or standard MR-AQCC.

Technically, after evaluating the ground state correlation energy (a single root AQCC calculation), the MR-AQCC-LRT procedure runs like a CI one: one obtains the eigenpair for all three excited states in one calculation, and there is no need for root following. No intruders appeared at least while obtaining three roots in both irreps.

In Table 1, the Q<sub>u</sub><sup>t</sup> values indicate that 2<sup>1</sup>A<sub>g</sub> and 2<sup>1</sup>B<sub>u</sub> states have covalent character, while the 1<sup>1</sup>B<sub>u</sub> state is ionic. The excitation energies obtained with AQCC and MR-AQCC-LRT for the covalent 2<sup>1</sup>A<sub>g</sub> and 2<sup>1</sup>B<sub>u</sub> states agree within 0.01 eV since the correlation energies for the ground state and these excited states agree within 0.1 eV. In contrast, the discrepancy is larger for the ionic state, with the AQCC-LRT excitation energy exceeding AQCC and MR-CISD+P by 0.4 eV. The underlying reason is that this state has a significantly higher correlation energy—nearly 2 eV—due to the absence of σ-correlation in the reference space. As a result, the AQCC-LRT correction is smaller, leading to a higher predicted excitation energy. Note that a significant improvement over standard MR-CISD is still observed. Moreover, for the 2<sup>1</sup>B<sub>u</sub> state the Q<sub>u</sub><sup>t</sup> value is reduced almost to half with MR-AQCC-LRT compared to MR-CISD, indicating reduced ionic character for this state. These results demonstrate that to avoid the unphysical mixing of covalent and ionic characters, a size-extensivity correction is essential in this case.

The transition densities of the relevant states are presented in Figure 4. As discussed previously,¹⁷,³⁹ the transition densities are located around the bonds for covalent states, whereas they are on the atoms for the ionic 1<sup>1</sup>B<sub>u</sub> state. More specifically, we observed the signature of σ-correlation,¹⁷,⁴⁴ as shown by additional transition density contributions near the atoms, which gives the individual atomic contributions the shape of d<sub>z<sup>2</sup></sub> orbitals. Such σ-correlation contributions can be visualized only by methods such as MR-CISD and MR-AQCC-LRT, both of which explicitly incorporate σ-correlation in the
<figure><img src="image_13.png" /><figcaption>**Figure 4.** Transition densities (isovalue: 0.002 au) between the ground and first three excited states of hexatriene, computed with the MR-AQCC-LRT method.</figcaption></figure>

wave function. The situation is different for second-order perturbation theory methods such as CASPT2. Such methods do provide an appropriate energy correction, but they usually do not give access to the perturbed transition densities that would be needed to visualize these σ-correlation contributions (and to provide appropriate perturbative corrections to computed transition properties).

The presented results on the intricate case of the three lowest excited states of all-trans hexatriene demonstrate that MR-AQCC-LRT is a viable tool for calculating transition properties, as the inclusion of the size-extensivity correction substantially improves the results compared to CI. Since the diagonalization of the same matrix yields the excited states, calculations are more straightforward than in the case of the original MR-AQCC. Additionally, analysis of the excited-state characters via transition densities is also possible. In the critical case of the ionic state, where the conditions for the LRT approach are not exactly met, it does not reproduce AQCC excitation energy precisely, but is still superior to MR-CISD.

**3.4. A First Step into GPU Application in COLUMBUS.**
Graphical Processing Units (GPUs) have significantly transformed the use of quantum chemical calculations.⁶¹ GPUs originated from applications in consumer video gaming and transitioned into systematic scientific programming with the advent of the Compute Unified Device Architecture (CUDA) framework.⁶² GPUs are specialized for highly parallel operations based on a Single-Instruction, Multiple-Thread (SIMT) architecture, in contrast to central processing units (CPUs), which process a sequence of threads with only a few tens of threads in parallel. The development of molecular quantum chemical GPU codes began with concentrating on the integral codes of GAMES-UK,⁶³ TeraChem,⁶¹ Gaussian 16,⁶⁴ PySCF,⁶⁵ to mention only a few programs for use in DFT. Taking TeraChem as an example, extensions to time-dependent density functional theory (TDDFT), CASSCF, and CASPT2 followed. Writing GPU-optimized programs can be a tedious task since two goals have to be achieved: (i) The algorithms must take advantage of the highly parallel GPU architecture, and (ii) the data transfer between CPU and GPU should be minimized. In TeraChem, the entire program is practically entirely based on GPU programming, where extensive sections of code had to be newly developed. Other codes followed a more conservative approach in transforming only selected segments of code to GPU usage.

The situation concerning different programming steps is more complex concerning the MR-CISD and the related MR-

AQCC methods in COLUMBUS. The sequence of individual steps is illustrated in a block diagram in Figure 5. The wave

<figure><img src="image_14.png" /><figcaption>**Figure 5.** List of individual steps for a single-point MR-CISD/AQCC calculation.</figcaption></figure>

function optimization step (determining the lowest eigenvalues and corresponding eigenvectors of the Hamiltonian matrix) is typically by far the most time-consuming. It is this step we are concentrating on in this work with the GPU implementation.

In short, the wave function optimization step in COLUMBUS is based on an uncontracted MR-CISD⁵⁰ calculation using GUGA⁷ with a direct CI technique.⁶⁶ All CSFs are represented in a distinct row table (DRT) shown graphically as a Shavitt graph⁷,⁶⁷ in Figure 6. This graph is divided into two parts, with the external orbitals at the bottom and the internal orbitals at the top. The external orbitals are limited to a cumulative occupation of up to two electrons, corresponding to the MR-CISD wave function expansion. The external part of the Shavitt graph is much simpler and regular, whereas the internal part can have any desired complexity. Once an internal walk is included in the expansion space, it is always combined with all possible external walks. This results in efficient computational kernels consisting of dense matrix–matrix products, matrix-vector products, and vector dot products with respect to the external orbitals,⁸,⁶⁸,⁶⁹ which may be offloaded onto accelerators and GPU hardware. It is noted that with typical modern basis sets, the number of external orbitals is typically much larger than the number of internal orbitals. The walks through the DRT (the CSFs) can be classified according to the passage through one of the four nodes, Z, Y, X, and W, at the interface between internal and external orbitals. Z-walks contain the reference CSFs and other configurations with zero excitation in the external space, Y-walks characterize single excitations, X-walks double excitations where the two external electrons are triplet coupled, and in W-walks the two external electrons are singlet coupled.

The dominant operation within the Davidson subspace method is the calculation of the matrix-vector product

$$w = Hx \qquad (16)$$

of the Hamiltonian matrix H with an expansion vector x of the Davidson procedure.⁷⁰,⁷¹ The matrix-vector products in eq 16 are computed from the underlying one- and two-electron Hamiltonian integrals ($h_{pq'} g_{pqrs}$) and the coupling coefficients ($\langle\tilde{m}'|E_{pq}|\tilde{m}\rangle$, $\langle\tilde{m}'|e_{pqrs}|\tilde{m}\rangle$). The $\tilde{m}$ denote the different walks (CSFs) through the DRT (Figure 6), and the $E_{pq}$ and $e_{pqrs}$ are one- and two-electron excitation operators. The general equation for this operation (eq 16) is
<figure><img src="image_15.png" /><figcaption>Figure 6. Graphical representation of a typical DRT (Shavitt graph) representing the CSF expansion space for a given set of electrons, spin, and number of orbitals. The lines in bold designate reference configurations. Each distinct, strictly ascending path from the origin to the top describes one of the CSFs of the expansion in terms of cumulative spin-coupling of the orbitals.</figcaption></figure>

$$w_{m',k} = \sum_{m,pq} h_{pq} \langle \tilde{m}' | E_{pq} | \tilde{m} \rangle x_m + \frac{1}{2} \sum_{m,pqrs} g_{pqrs} \langle \tilde{m}' | e_{pqrs} | \tilde{m} \rangle x_m \quad (17)$$

The coupling coefficients are computed on the fly as necessary. These coupling coefficients are very sparse, so the procedure adopted in the MRCI code involves the efficient computation of only the nonzero coupling coefficients. The coupling coefficients are computed in blocks based on the partitioning of orbitals and the types of CSFs, W, X, Y, or Z, associated with the respective Hamiltonian block type. The whole computation of eq 17 is driven by the internal indices and the coupling coefficients belonging to a given internal index set. Only the contribution of internal orbitals is computed in this step. The contribution of the external orbitals is calculated on the fly, leading to the matrix multiplication operations shown in Figure 7 (see also refs 8,68). For a demonstration of the GPU implementation, we focus on this equation in the most critical section, which is given by walks through the interface nodes X or W leading to the XX, WW, and WX pairs.

The modular treatment of matrix multiplication kernels via standard subroutine calls to the Basic Linear Algebra Subprograms (BLAS) library was implemented in COLUMBUS some time ago⁸ to facilitate the adaptation of the COLUMBUS code to the evolution of hardware and software

w = Hx Update Scheme

<figure><img src="image_16.png" /><figcaption>Figure 7. Update scheme for the vector w (eq 16) for the special case of CSFs with two external indices and a representative choice of two-electron integral contributions as an example. Operations involving the GPU are in red.</figcaption></figure>

components. This strategy is also beneficial for GPU programming. As will be described below, only a few changes in the source code need to be made, such as exchanging the CPU-based BLAS library with a corresponding GPU-optimized CUDA library (cuBLAS) and inserting commands taking care of the (limited) data transfer to and from the GPU.

Figure 7 illustrates this strategy for one representative example. It is shown that with one upload of w and x vector segments and the transfer of one block of integrals belonging to a fixed internal index pair kl, a large number of matrix multiplication updates of the vector w (eq 16) can be made coupling all CSFs according to the GUGA formalism. These operations can be performed independently and asynchronously without any further data transfer between CPU and GPU. Thus, only a relatively small effort is required for data transfer of two-electron integrals and CI vector segments, as a large number of matrix multiplications can be performed on the GPU using these data. The number of these multiplications will increase strongly with the increase of the active space since an increasing number of coupling elements will be produced for a set of integral indices.

The GPU implementation is currently based on CUDA-Fortran as provided in NVIDIA's Compiler suite (NVIDIA HPC SDK).⁶² Variables (matrices and vectors) residing on the GPU are allocated directly using Fortran declarations and variable-to-variable copying, avoiding the need for explicit CUDA function calls. In the current implementation, only device initialization (creating the cuBLAS handle) and GPU call synchronization are performed using dedicated CUDA function calls. These calls can be easily replaced in the future with corresponding functions from other implementations, allowing for straightforward porting of the code to other hardware accelerators that use different programming models. The matrix multiplication itself is performed by calling the corresponding vendor-provided Linear Algebra PACKage (LAPACK) DGEMM routine through a wrapper function. The overhead caused by calling the corresponding BLAS routines on the GPU will be compensated by (a) the speedup of the DGEMM routine for large matrices on the GPU compared to the CPU implementation and (b) the fact that GPU calls can be run asynchronously, meaning the CPU can continue processing the loop over all internal contributions while the multiplication on the GPU is still in progress.

An implementation for handling the two external integrals, as shown in Figure 7, has been recently performed, and
<figure><img src="image_17.png" /><figcaption>**Figure 8.** Diagram showing the implementation of projected openCAP in COLUMBUS. Initially, a CAP-free MRCI calculation is performed to obtain a limited number of eigenstates. The one-electron densities, $γ_{II}$ and transition densities, $γ_{IJ}$ are used to calculate the CAP ($\hat{W}$) matrix elements. The complex Hamiltonian $\hat{H}(\eta)$ is diagonalized for several values of $η$ and the $η$-trajectory is used to obtain the Siegert energies.</figcaption></figure>

benchmarking tests are in progress. Treatment of the other integral types is straightforward and will follow analogous procedures. The current implementation is tightly integrated with COLUMBUS' approach to parallelization, utilizing the Global Arrays package in conjunction with an MPI-based interprocess communication. This enables the use of multiple GPUs in parallel, thereby mitigating the issue of a fixed, vendor-dependent memory limit standard to contemporary GPUs. The implementation was performed on the Sophia supercomputer at the Argonne Leadership Computing Facility. Preliminary versions of the implementation can be found in a dedicated branch of COLUMBUS' GitLab repository.

**3.5. Complex Absorbing Potentials for Electronic Resonances.** The ability to describe autoionizing states accurately remains a challenging problem in the electronic structure community. Electronic resonances, a type of autoionizing state, are formed through the resonant attachment of electrons in low-energy electron (LEE)-induced reactions. LEE-induced reactions are prevalent in many areas, including astrochemistry, radiation chemistry, solid-state processes, and others.⁷²⁻⁷⁶

Electronic resonances are categorized into either single-particle (1p) or two-particle one-hole (2p-1h) resonances. A single-particle shape resonance is formed when a centrifugal potential holds the excess electron bound to the neutral system. These resonances decay through a one-electron process and are characterized by shorter lifetimes.⁷⁷ In 2p-1h resonances, the excess electron is bound to an excited state of the neutral target and, based on the decay mechanism, can be further categorized into “core-excited shape” and Feshbach.⁷⁷ Most studies on electronic resonances in the literature have been limited to shape resonances. Fewer studies focus on Feshbach resonances despite their presence in a wide variety of electron-driven reactions.⁷⁸⁻⁸⁰ Additionally, resonant-channel coupling between shape and Feshbach channels facilitates electron-driven processes in nucleobases and similar aromatic systems.⁸⁰⁻⁸² A multireference treatment of Feshbach resonances and their mixing with shape channels is pertinent for their accurate characterization.⁸¹,⁸³ To address this deficiency, an approach based on complex absorbing potentials was recently implemented by creating an interface between

COLUMBUS and the OpenCAP software.⁸⁴ COLUMBUS offers several advantages, and its flexibility in describing the multireference space is crucial for describing resonances involving multireference character. Furthermore, the availability of analytic gradients and nonadiabatic couplings makes it an attractive package for future developments involving methods for resonances. The basics of this development are summarized here.

While there are many approaches to describe resonances, we have focused on non-Hermitian quantum mechanics⁸⁵ techniques, where resonances can be obtained as complex solutions to the non-Hermitian eigenvalue problem, with complex (or Gamow–Siegert) eigenvalues of the form

$$E_{\text{res}} = E_{\text{R}} - i\Gamma/2 \qquad (18)$$

where $E_R$ is the energy or position of the resonance state, and $\Gamma$ is its associated width. Inverse width provides information on the lifetime of the resonance state. In CAP, a non-Hermitian Hamiltonian is obtained by adding a complex local potential to the physical Hamiltonian for a given system⁸⁶

$$\hat{H}(\eta) = \hat{H} - i\eta\hat{W} \qquad (19)$$

where $\hat{W}$ is a one-particle potential, and $\eta$ is the CAP strength. A commonly employed one-particle potential in CAP calculations is a quadratic potential (also otherwise known as box-potential) of the form⁸⁶

$$\hat{W}_k = \begin{cases} 0, & r_k < r_k^0 \\ (|r_k| - r_k^0)^2, & r_k > r_k^0 \end{cases} \qquad (20)$$

where $k = x, y, z$ and $\hat{W} = \sum_k \hat{W}_k$. The CAP onset makes sure that the interaction region of the resonance wave function remains unperturbed from the presence of the one-particle potential.

A critical aspect of CAP is the strong parametric dependence of the complex energies on the CAP strength ($\eta$). The real and complex part of the eigenvalues as a function of $\eta$ ($\eta$-trajectory) tend to spread out in the complex plane, and the stabilization point is revealed by minimizing the first-order term in the Taylor series expansion of the $\eta$-dependent
energies. The complex energies corresponding to the stabilized point identify the optimum value.⁸⁶

$$|E(\eta) - E(0)| = \left| \eta \frac{dE}{d\eta} \right| + O(\eta^2) \qquad (21)$$

The sensitivity of this parametric dependence on η can be further minimized. Many ways have been prescribed for the secondary correction,⁸⁶,⁸⁷ and one of the straightforward ways is to minimize the second-order correction in the Taylor expansion as proposed by Riss and Meyer⁸⁶

$$|E(\eta) - E(0) - \eta \frac{dE}{d\eta}| = \left| \eta^2 \frac{d^2 E}{d\eta^2} \right| + O(\eta^3) \qquad (22)$$

CAP can be introduced into the electronic structure calculation at any stage. The conventional route involves adding the CAP at either the Hartree–Fock (HF) step or the correlated step, depending on the given CAP strength. These calculations are computationally expensive and require modifications to the traditional framework to accommodate changes resulting from the complex algebra involved. Additionally, determining the optimum resonance parameters requires complex-valued energies for multiple values of the CAP strength parameter. Therefore, multiple electronic structure calculations must be carried out to analyze the η-trajectories. An alternative to conventional CAP is the projected CAP, where CAP is projected onto a correlated subspace.⁸⁸,⁸⁹

We have implemented projected CAP with the MRCI code in COLUMBUS, as outlined in Figure 8. The OpenCAP program by Gayvert and Bravaya was used for CAP calculations.⁸⁴ Details on the projected CAP-MRCI implementation can be found in refs 90,91. This implementation was tested on several systems and successfully produced complex-valued potential energy surfaces, as well as demonstrated mixing between shape and Feshbach resonances. This approach demonstrated the profound effect the mixing has on the lifetimes of the resonances, with potential implications in the behavior of nucleobases under LEE damage and other relevant applications.⁹⁰,⁹¹

**3.6. Hybrid Antisymmetrized Coupled Channels.** A customized code version of COLUMBUS delivers the CI functions for applications in attosecond photoemission from molecules in strong laser fields.⁹²,⁹³ That process is highly nonperturbative and imposes special requirements on the representation of the system. Emission rates and spectra are susceptible to thresholds and initial and final state correlations.⁹⁴ Therefore, correlated ground and excited states must be accurately represented. In addition, photoelectrons are emitted into a vast spectrum of continuum states, from threshold up to hundreds of eV. For single-electron emission, the CI space is augmented by a complete single-electron basis represented on a numerical grid. An ansatz is made for the time-dependent wave function as a linear combination of one or more neutral states and “channel” functions of ionic states that are augmented by the numerical basis

$$|\Psi, t\rangle = \sum_N |N\rangle c_N(t) + \sum_A \left\{ \sum_i a_i^\dagger |A\rangle c_i^A(t) + \sum_\alpha a_\alpha^\dagger |A\rangle c_\alpha^A(t) \right\} \qquad (23)$$

with complex time-dependent coefficients $c_N(t)$, $c_i^A(t)$ and $c_\alpha^A(t)$. In this “hybrid” representation, the neutral states $|N\rangle$ and the ionic states $|A\rangle$ are generated by COLUMBUS in MCSCF-MRCI, based on the molecular orbitals $|li\rangle$ of the neutral state at its equilibrium geometry. Using the orbitals at the neutral positions also for the ions is justified in the approximation that relaxation of the nuclear positions occurs on a longer time scale than the electronic emission. The creation operators $a_i^\dagger$ generate all molecular orbitals. The $a_\alpha^\dagger$ generate numerical basis functions $\{|\alpha_\perp\rangle\}$ that are orthogonalized to the orbital space $\{|li\rangle\}$ as

$$|\alpha_\perp\rangle = \left( 1 - \sum_i |li\rangle\langle li| \right) |\alpha\rangle \qquad (24)$$

When calculating the effect of two-particle operators $\hat{T} = T_{k,l,mn}a_k^\dagger a_l^\dagger a_m a_n$ in this hybrid basis, one needs up to three-particle density matrices of the ions w.r.t. the molecular basis

$$\langle a_i^\dagger A a_k^\dagger a_l^\dagger a_m a_n a_j^\dagger B \rangle = \langle A a_i a_k^\dagger a_l^\dagger a_m a_n a_j^\dagger B \rangle \qquad (25)$$

and similarly Dyson-orbital-like objects with up to 5 indices for matrix elements between the $|a_i^\dagger A\rangle$ and $|N\rangle$. At present, these are generated by using an (unpublished) version of COLUMBUS that outputs the $|N\rangle$ and $|A\rangle$ as lists of determinants and coefficients. The matrix elements of single- and two-particle operators involving the Gaussian-based orbitals $|li\rangle$ and the numerical basis $|\alpha\rangle$ are computed by numerical quadrature. Linear dependence of the basis is controlled by a low-rank representation of the corresponding pseudoinverse using a Woodbury-type formula.

The time-dependent Schrödinger equation is solved with the ansatz (23) with absorbing boundary conditions. Photoemission spectra are extracted by analyzing outgoing flux using the t-Surff method,⁹⁵ which accounts for the fact that emission spectra are modified by the time-dependent field outside the simulation volume.

Figure 9 shows strong-field photoemission spectra of CO₂ in an elliptically polarized field at a wavelength of 800 nm (first published in ref 96). Neutral and up to 6 ionic states were used in this calculation. The strong field dominates, and the spectrum does not exhibit any multipolar characteristics. The process is in the intermediate region between tunneling in the strong field and extreme multiphoton emission. Tunneling

<figure><img src="image_18.png" /><figcaption>**Figure 9.** Strong field photoemission from a CO₂ molecule by a laser pulse with wavelength 800 nm and elliptic polarization. The molecule is aligned along the Z-axis, ellipticity is $\epsilon = 0.87$. The polarization plane is XY in (a) and XZ in (b–d). Energies extend to 1.5 au, and the cut plane is indicated in the subfigures.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 10.** Comparison between MR-AQCC $N_U$ values and FT-DFT $N_{FOD}$ numbers with M05-2X density functional for structures with singlet ground electronic state (structures 1–16, 18, 25–29 of ref 110) using (a) the original $T_{el}$ of 16,200 K and (b) the improved $T_{el}$ of eq 31 of 12,200 K. Reproduced with permission from ref 110. Copyright 2023 the Royal Society of Chemistry.</figcaption></figure>

produces energetically broad distributions modulated by ripples due to multiphoton processes. The effects of ellipticity are most visible in the cuts through the polarization plane in Figure 9c. Symmetries of the HOMO are reflected in corresponding nodes of the photoemission spectra.

## 4. APPLICATIONS: OPEN-SHELL SYSTEMS

### 4.1. Calibration of FOD Parameters Using MR-AQCC Unpaired Densities.
Graphene and graphene-based nano-materials form the basis for numerous interesting applications, including energy storage devices, optoelectronics, flexible conductive films for displays and electrodes, and anticorrosion technology.⁹⁷,⁹⁸ When broken down into nanoscale fragments, PAHs⁹⁹ are formed, exhibiting notable and tunable diradical and polyradical characteristics.¹¹,¹⁰⁰⁻¹⁰³ These properties make PAHs attractive for spintronics, optoelectronics, conductivity, and magnetism,¹⁰⁴,¹⁰⁵ though their complex electronic structure poses significant characterization challenges. The characterization of these molecules, particularly in terms of their biradical character, provides valuable information for developing new PAH-based materials. The number of effectively unpaired electrons ($N_U$),¹⁰⁶⁻¹⁰⁸ accessible through multireference (MR) methods,²⁰,⁵⁰ provides valuable insights into polyradical systems. An alternative descriptor based on DFT calculations is the fractional occupation number-weighted electron density (FOD)⁶,¹⁰⁹ obtained from finite-temperature DFT (FT-DFT) without fictitious temperature-dependent energy functionals. This method effectively simulates static correlation straightforwardly based on the electronic temperature ($T_{el}$), which acts as a parameter. Although multireference methods can provide reliable information about the open-shell character of PAHs, they are relatively expensive, and DFT-based methods are preferred for large-scale, routine applications. However, in the case of FT-DFT applications, the electronic temperature parameter $T_{el}$ needs to be adjusted. In the present report, this adjustment is summarized for a series of functionals through fitting MR-AQCC $N_U$ values.¹¹⁰,¹¹¹

For optimizing the $T_{el}$ parameter, 22 PAHs ranging from acenes, indenoacenes, zethrenes, and other biradicaloid PAHs¹¹⁰ were selected and benchmarked using MR-AQCC calculations. The $N_U$ values and the corresponding unpaired electron densities¹⁰⁶,¹¹² were obtained using Head-Gordon's nonlinear formula¹⁰⁷

$$N_U = \sum_{i=1}^{N} n_i^2 (2 - n_i)^2 \qquad (26)$$

Here, N is the number of natural orbitals (NOs), and $n_i$ is the occupation of the ith NO. This formula emphasizes orbitals with occupations close to one while suppressing contributions near 0 and 2. The geometries of all 22 structures were computed using the TPSS¹¹³ functional and the def2-TZVP basis set¹¹⁴ and are available in the Supporting Information of ref 110.

In the FOD analysis,⁶ the fractional orbital density is defined as

$$\rho^{\text{FOD}} = \sum_{i=1}^{N} (\delta_{1} - \delta_{2}f_{i})|\phi_{i}(r)|^{2} \qquad (27)$$

where $\delta_1$ and $\delta_2$ are constants. These are both set to unity if the orbital energy is below the Fermi energy ($E_F$). $\delta_1$ is set to 0 and $\delta_2$ to -1, respectively, if the orbital energy is above $E_F$. The $\phi_i$ terms represent molecular spin orbitals, and the summation extends over all orbitals N. The $f_i$ values are derived from the Fermi-Dirac distribution, providing fractional orbital occupation numbers as

$$f_i = \frac{1}{e^{(\epsilon_i - E_F)/kT_{el}} + 1} \qquad (28)$$

where $\epsilon_i$ represents the orbital energy. The number of “hot” electrons, $N_{FOD}$, is defined as

$$N_{\text{FOD}} = \sum_{i=1}^{N} (\delta_{1} - \delta_{2}f_{i}) \qquad (29)$$

In the original work,⁶ the electronic temperature ($T_{el}$) was determined as a function of the percentage of nonlocal Fock exchange admixture ($a_x$) in the functional

$$T_{el} = 20000 \text{ K} \times a_x + 5000 \text{ K} \qquad (30)$$

The plot in Figure 10a shows that for the M05-2X functional ($a_x$ = 0.56), there is a good linear correlation between $N_{FOD}$ values and $N_U$ values, using the original $T_{el}$ parametrization of eq 30. However, the original $T_{el}$ value is much too large leading to a strong overestimation of the $N_{FOD}$ values in comparison to respective $N_U$ data. Optimization of $T_{el}$ in steps of 1000 K and refinement in steps of 200 K, the
<figure><img src="image_21.png" /><figcaption>Figure 11. Comparison between MR-AQCC $N_U$ values and FT-DFT $N_{FOD}$ numbers for the ωB97XD density functional for structures with singlet ground electronic state (structures 1–16, 18, 25–29 numbered as in ref 111) using (a) $T_{el}$ of eq 31 of 8530 K and (b) the improved $T_{el}$ of eq 32 of 14,000 K. Reproduced with permission from ref 111 under a Creative Commons Attribution 4.0 International License. Copyright 2024 Carvalho et al.¹¹¹</figcaption></figure>

slope is improved, and at 12,200 K, a slope close to one is achieved (Figure 10b).

Similar investigations with the TPSS, M06–2X, and B3LYP functionals revealed that better agreement could be achieved using the improved equation for these classes of functionals as well¹¹⁰,¹¹¹

$$T_{el} = 10762 \text{ K} \times a_x + 6140 \text{ K} \quad (31)$$

Further tests indicated that this improvement extends to the double-hybrid functional B2PLYP as well.¹¹¹

Range-separated hybrid functionals have also been investigated.¹¹¹ Application of eq 31 developed for GGA functionals to the range-separated ωB97XD functional ($a_x = 0.222$) leads to a $T_{el}$ of 8530 K. Although the correlation between the corresponding $N_{FOD}$ and $N_U$ values shows good linearity (Figure 11a), the slope is much too small. By increasing $T_{el}$ to 14,000 K, an excellent overall correlation is achieved (Figure 11b).

The same adjustments of $T_{el}$ were performed for the range-separated functionals ωB97M-V, CAM-B3LYP, LC-ωPBE, and MN12-SX. It was found that the resulting $T_{el}$ values could be fitted against the $β$ parameter, which is the difference between the HF exchange in the long-range limit ($α + β$) and the corresponding percentage in the short-range limit¹¹⁵ ($α$) leading to the expression.¹¹¹

$$T_{el} = 6650 \text{ K} \times β + 9010 \text{ K} \quad (32)$$

In addition to the parametrization of $T_{el}$ for various types of DFT functionals, an optimal Fermi temperature has also been determined for the semiempirical GFN2-xTB method,¹¹⁶,¹¹⁷ allowing reliable $N_{FOD}$ data at even much faster computational speeds.¹¹¹ In summary, a set of Fermi temperature parameters $T_{el}$ has been derived for a range of DFT functionals based on a larger selection of open shell PAH compounds, leading to $N_{FOD}$ values that reproduce the MR-AQCC $N_U$ values very well. This opens up excellent opportunities for computationally efficient predictions of the biradical character, especially for larger molecules. So far, the investigations have focused on PAHs with applications also to B,N doped periacenes (see also Section 5.2) and CC dissociation processes in ethane, ethylene and acetylene.¹¹¹ In all these cases good agreement with the number of unpaired electrons $N_U$ were found. Thus, we expect that carbon biradicals and similar systems should be well covered by the present $T_{el}$ parametrization. To demonstrate its

applicability beyond this set of compounds, more extended classes of molecular structures still need to be investigated.

## 4.2. Electronic States in Single-Vacancy (SV) Graphene Defects.
Graphene is renowned for its exceptional electronic, thermal, and mechanical properties.⁹⁷ Doping of graphene sheets is frequently used as an efficient tool to modulate their electronic properties. One important class of graphene defects is carbon vacancies, which can occur either naturally due to growth defects or by collisions with high-energy species.¹¹⁸ The removal of a carbon atom from the regular graphene network creates dangling bond carbons, leading to low-lying electronic states.¹¹⁹,¹²⁰ The detailed structural characterization of these states and their interplay with geometrical changes comprise interesting but challenging questions. Reliable approaches require multireference methods (MR) that describe both the ground and excited states equally well. MR-CISD,¹²,⁵⁰,⁵¹ as available in COLUMBUS, provides a good basis for such investigations. This contribution summarizes our previous work on a single vacancy defect in pristine pyrene.¹¹⁹⁻¹²² By removing one carbon atom from the center of pyrene (Figure 12a) leading to pyrene-1C, three

<figure><img src="image_22.png" /><figcaption>Figure 12. (a) Pristine pyrene, (b) unrelaxed pyrene-1C, carbon dangling-bond centers are shown in black, (c) pyrene-1C after geometry relaxation.</figcaption></figure>

dangling bonds are formed upon rigid C removal (Figure 12b). On planar geometry relaxation, two interact to form a covalent bond, leaving only one dangling carbon bond (Figure 12c). Planar pyrene-1C has a quasi-degenerate ground state with three low-lying excited states. This finding prompted us to investigate further whether out-of-plane distortions could result in additional stabilization and how excited states could be affected.¹²²
The geometry of pyrene-1C was fully optimized in C₂ᵥ symmetry using a CASSCF(16,16) approach that included all π-orbitals and the carbon dangling bond for each of the four lowest states (¹A₂, ³A₂, ¹B₁, and ³B₁). The CAS was then systematically reduced to a CAS(6,6) and finally to a compact CAS(4,3) that was able to reproduce the CAS(16,16) excitation energies and geometries well. Final calculations were performed at SA4-CASSCF(4,3) (averaged over the four lowest states) and MR-CISD+Q(4,3) (+Q indicates Davidson correction¹⁵,⁵¹) levels using the 6–31G* basis set.¹²³ All stationary structures described in this work were characterized by harmonic frequencies calculated at the SA2-CASSCF(4,3)/6–31G* level.

The C₂ᵥ structures optimized for the ¹A₂, ³A₂, ¹B₁, and ³B₁ states were all found to be transition states, with ³B₁ being more stable than the others by less than 0.2 eV. Following the imaginary modes, four nonplanar optimized structures with C₅ symmetry were obtained, of which states 1A' and 3A'' were transition state structures, and ¹A'' and ³A' were local minima. Further distortion of ¹A' from C₅ to C₁ symmetry resulted in a ¹A structure, whereas the distortion of ³A'' led to a ³A' structure. An energy diagram of these structures is presented in Figure 13. Among these structures, the L-shaped conformation

<figure><img src="image_23.png" /><figcaption>Figure 13. Energies ΔE (eV) are relative to the ³B₁ state using the MR-CISD+Q(4,3) method for the optimized geometries of the ¹A₂, ³A₂, ¹B₁, and ³B₁ states in C₂ᵥ symmetry, the optimized geometries of the ¹A', ³A', ¹A'', and ³A'' states in C₅ symmetry and for the ¹A and ³A states in C₁ symmetry. Transition states are labeled in red.</figcaption></figure>

of the ¹A' state, obtained from the distortion of the ¹B₁ state, was considerably more stable than the other distorted C₂ᵥ structures. Its stabilization amounted to 1.13 eV in comparison to the ³B₁ structure. The higher stability of this conformation is derived from a 3c-2e bond formed between the carbon of the dangling bond and a C–C bond of the cyclopentadienyl fragment. Despite the stability of this structure, it was still not a local minimum on the potential energy surface. Further reduction of the symmetry to C₁ allowed the localization of the 3c-2e bond into a σ-bond between the dangling of the carbon bond and a carbon of the cyclopentadienyl ring. This structural rearrangement led to a highly distorted structure, which, nevertheless, is about 3.77 eV lower in energy compared to the reference ³B₁ structure because of the

formation of a CC σ-bond. Calculating a triplet state using the ¹A geometry as a starting point resulted in a similar geometry but 1.55 eV higher in energy.

Vertically excited states for the symmetries B₁ and A₂ of the lowest singlet states are shown in Figure 14a, computed at the C₂ᵥ geometry of the ³B₁ structure. There are three excited states below 1 eV, while the remaining are all above 4 eV relative to the ground state. The most intense excitation, calculated from the ¹B₁ state, is to the 4¹A₂ state, with an excitation energy of about 7.8 eV. A similar pattern is observed for the other C₂ᵥ and C₅ optimized structures. The only exception is the ¹A' optimized structure, which has a closed-shell configuration due to the formation of a 3c-2e bond, resulting in a wide gap between the ground state and the first three excited states, which are about 3 eV higher in energy (Figure 14b). The remaining excited states are all above 5 eV. Two intense excitations are observed from the 1 ¹A' state: one to 3¹A'' and the other to 4¹A''.

A systematic distortion of the planar pyrene-1C structure resulted in several stationary points from which four structures were minima. A highly distorted closed-shell ¹A structure with C₁ symmetry was the most stable one, with 3.77 eV more stable than the planar ³B₁ reference structure, followed by an open-shell ³A structure, which has a similar geometry to ¹A but is 1.55 eV higher in energy. The other minima (¹A'' and ³A'), which are less distorted compared to the planar structure but have a dangling bond, are at most 0.02 eV more stable than the original ³B₁ reference structure. In larger sheets, greater restraining forces will exist, certainly disfavoring the C₁ symmetry structure. The others should be more adaptable in their current structures. Investigations into larger graphene sheets will be necessary to resolve this question conclusively.

### 4.3. Construction of Polyradical Character with Phenalenyl Oligomers.
Carbon-based radicals have recently been the subject of intense study due to their unique electronic and magnetic properties, as well as their potential applications in materials chemistry and biosciences.¹²⁴ The high chemical reactivity of radicals presents a considerable challenge for their synthesis and characterization, often necessitating a combination of solution and on-surface chemistry.¹²⁵ Among these, the structures composed of benzenoid rings arranged in triangular patterns are regarded as promising materials in spintronics applications,¹²⁶⁻¹²⁸ and for the design of thermally active delayed emitters.¹²⁹

A small singlet–triplet (S₀/T₁) energy gap (ΔE) is the main prerequisite for a desirable magnetic switching in these materials.¹³⁰ Despite their reactivity originating from open-shell electronic structures, triangle-shaped systems, and their two-dimensional extensions have been synthesized and characterized using advanced on-surface approaches.¹²⁸,¹³¹⁻¹³⁴ Studies of triangulene dimers, directly coupled or linked via bridging groups such as acetylenes and phenyl, demonstrate that the nature of these connections can tune their magnetic properties and ground-state spin states well.¹²⁷,¹³⁴⁻¹³⁹ The triangle-shaped radical compounds and their extensions form alternant π-conjugated systems and belong to the class of non-Kekulé structures. In these structures, the spin quantum number (S) of the ground state is predicted according to Ovchinnikov’s rule¹⁴⁰ and Lieb’s theorem¹⁴¹ as

$$S = \frac{N^* - N}{2} \qquad (33)$$
<figure><img src="image_24.png" /><figcaption>**Figure 14.** Vertical excitation energies $\Delta E$ (eV) relative to (a) the optimized $^3B_1$ structure and (b) to the $^1A'$ structure computed at the MR-CISD +Q(6,6) level. Oscillator strengths refer to the singlet states (represented by the red vertical lines). For triplet states, black vertical lines only indicate the energetic position of the electronic state. Adapted with permission from ref 122. Copyright 2023 American Chemical Society.</figcaption></figure>

where $N^*$ and $N$ are the number of starred and unstarred sites in the alternant $\pi$-conjugated systems, respectively.

In this contribution, we report on an initial study on covalently linked dimers of phenalenyl, the smallest member of triangle-shaped molecules with highly delocalized radical structure (Figure 15a), high electrical conductivity, and

<figure><img src="image_25.png" /><figcaption>**Figure 15.** (a) Phenalenyl: $\alpha$ and symmetry equivalent C atoms denote positions with SOMO contributions, whereas $\beta$ positions do not contribute to the SOMO, (b) phenalenyl dimers **1–3** linked according to $\alpha/\beta$ positions; spin multiplicity given according to Ovchinnikov's rule.</figcaption></figure>

distinctive magnetism.<sup>142–146</sup> The $\alpha$ positions label those atoms where the spin density in the singly occupied molecular orbital (SOMO) is located, whereas $\beta$ atoms do not carry spin density. Due to their size, phenalenyl dimers are ideal for studying how the different linking types between phenalenyl units influence their electronic properties. Depending on the coupling sites, the following structural isomers can be identified: (i) isomers coupled with the spin positions ($\alpha,\alpha$-dimer **1**, Kekulé structure), nonspin positions ($\beta,\beta$-dimer **2**, non-Kekulé structure), and spin-nonspin positions ($\alpha,\beta$-dimer **3**, non-Kekulé structure), illustrated in Figure 15b. All three structures are alternant $\pi$-conjugated systems. Applying Ovchinnikov's rule to these systems suggests a singlet ground state for $\alpha,\alpha$- and $\beta,\beta$-dimers, while the ground state of the $\alpha,\beta$-dimer should be a triplet spin state. This contribution discusses the electronic structure of the three isomers linked via an acetylene bridge. For this purpose, extended multireference calculations were performed to evaluate the singlet–triplet splitting energy in comparison with the prediction by Ovchinnikov's rule and to characterize the open shell character of the dimers employing numbers of unpaired electrons ($N_U$).

The calculations were performed using the MR-AQCC<sup>20</sup> method in combination with the orbitals obtained from a state-averaged (SA)-complete active space self-consistent field theory consisting of eight electrons in eight orbitals (SA-CASSCF(8,8)). The same CAS was used as reference space in the MR-AQCC calculation. The reference space has been chosen based on the NO occupations by including the main open-shell orbitals within an occupation range of ~1.90 e to ~0.1 e. The $\sigma$ space orbitals were kept frozen after the SCF level, while all $\pi$ orbitals were included, as has been successfully used for a series of PAHs in previous investigations.<sup>11,100,147</sup> All calculations were performed with the 6–31G(d) basis set.<sup>148</sup> The singlet structures were optimized at the unrestricted density functional theory (UDFT) level using the B97-D functional<sup>149</sup> and the TZVP basis set.<sup>114</sup> S-T splittings were calculated as $\Delta E_{ST} = E_T - E_S$, computed as the vertical energy difference between the lowest singlet ($E_S$) and triplet ($E_T$) states, where the positive values indicate that the singlet is the lower energy state. A measure for the total number of unpaired electrons ($N_U$) was calculated using eq 26.

In agreement with the results on directly coupled dimers<sup>138</sup> and according to Ovchinnikov's rule (eq 33) and Figure 15b, the number of star and nonstar sites assigns a singlet spin state to dimers **1** and **2** and a triplet spin state to dimer **3**, assuming that the acetylene bridge is involved in conjugation. S-T splitting energies are collected in Table 2. A large positive S-T

<table><caption>**Table 2.** Singlet-Triplet Splitting Energies and $N_{U,S}$ (Singlet State) and $N_{U,T}$ (Triplet State) Values for 1-3 Computed at MR-AQCC Level</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>connection</td><td>α,α</td><td>β,β</td><td>α,β</td></tr><tr><td>ΔE<sub>ST</sub> (eV)</td><td>0.924</td><td>-0.003</td><td>-0.004</td></tr><tr><td>N<sub>U,S</sub> (e)</td><td>1.66</td><td>2.70</td><td>2.61</td></tr><tr><td>N<sub>U,T</sub> (e)</td><td>2.58</td><td>2.70</td><td>2.63</td></tr></tbody></table>

splitting value is found for **1** ($\Delta E_{ST} = 0.92$ eV). It is related to a strong interaction due to the involvement of $\alpha$ atoms as connectors between the two phenalenyl rings. This case aligns well with Ovchinnikov's rule, which predicts that **1** should have a singlet ground state. A slightly negative value is observed for **3** ($\Delta E_{ST} = -0.004$ eV), predicting a quasi-degeneracy of the singlet and triplet states. Formally, the prediction by the Ovchinnikov rule is correct but misses the quasi-degeneracy. A tiny negative splitting is found for **2** ($\Delta E_{ST} = -0.003$ eV),
<figure><img src="image_26.png" /><figcaption>Figure 16. Unpaired electron densities (isovalence 0.004 e/Å³) for the singlet and triplet states of 1–3.</figcaption></figure>

which implies a quasi-degeneracy as well as slightly favoring the triplet state again. This contradicts the prediction based on the Ovchinnikov rule. The $N_{U,S}$ values (Table 2) are significantly larger for both **2** and **3** than for **1**. They are about twice the $N_U$ value of phenalenyl (1.36 e) and indicate weak interactions between the phenalenyl units in these two cases. $N_{U,T}$ values are almost identical to the singlet ones for structures **2** and **3**. The stronger interaction between the two phenalenyls in **1** reduces the $N_{U,T}$ value somewhat.

Figure 16 displays the calculated unpaired densities for the three isomers, showing significant differences in the distribution of unpaired electrons in the acetylene bridge. The unpaired density for the singlet state is intensely concentrated in the $(\alpha,\alpha)$-type bridge in **1** (Figure 16a), reducing the radical character of the phenalenyl rings. On the contrary, the other two connection types $(\beta,\beta)$ and $(\alpha,\beta)$ (Figure 16b,c) do not lead to any (in **2**) or minimal (in **3**) unpaired densities in the acetylene bridge. The unpaired density for the triplet state of **1** (Figure 16d) shows a significant open-shell character localized in the phenalenyl rings, in contrast to the respective singlet case (Figure 16a). There is also substantial unpaired density localized in the acetylene bridge. The unpaired densities for **2** and **3** (Figure 16e,f) are practically identical to the open-shell singlet counterparts.

The results for different isomers of phenalenyl dimers connected by an acetylene bridge reveal the variability in their electronic structure and, consequently, their magnetic properties, depending on the type of connection. The connecting bridge plays a crucial role in these variations. Furthermore, it is suggested that the different roles of the bridge, particularly its interference with the radical character, may cause deviations from the spin predictions based on the Ovchinnikov rule in cases where the interactions between the two radical subsystems are weak due to the involvement of $\beta$-type connectors. Ovchinnikov's rule will always favor one multiplicity. In contrast, our results show that in some cases, singlet and triplet states can be quasi-degenerate, and both should be considered on equal footing. This situation highlights the importance of high-level computational studies performed at the multireference level.

### 4.4. Aromatic and Heteroaromatic Polyradicals.

Polyradicals emanating from aromatic and heteroaromatic rings are an interesting class of molecules; highly reactive due to the presence of dangling, unpaired electrons yet thermodynamically stabilized by the aromaticity present in the ring.¹⁵⁰ Such polyradicals can be formed via the electrocyclization of acyclic precursors¹⁵¹ or found in high-temperature combustion of cigarette smoke or alternative fuels.¹⁵² They are also vital intermediates in the fabrication of

polymeric materials.¹⁵³ The reactivity and fleeting existence of these species complicate their experimental characterization, while the complexity of their electronic structure necessitates careful theoretical analysis. For instance, using polyradicals containing two unpaired electrons as an example, two electrons in two orbitals give rise to six electron configurations. These electron configurations combine to form four electronic spin states: a triplet and three singlet states.³⁷,¹⁵⁴ Each of the singlet states requires more than one electron configuration (as represented by a Slater determinant) to fully capture the physical nature of each singlet state.

The proper characterization of these multiconfigurational singlet states requires multireference (MR) methods involving initial guess wave functions that contain more than one determinant. The multireference and complete active space (CAS) methods available in COLUMBUS are well-suited to describing multiconfigurational systems. Additionally, the energy manifolds of polyradicals contain a high density of closely lying states. In these instances, the state-averaging available in COLUMBUS enables the development of a state-averaged, multireference wave function that is optimized in an averaged manner for all states. That reference wave function is ideal for subsequent single-state analysis using the highly correlated levels of theory available in COLUMBUS, such as MR-CISD and MR-AQCC. The parallel implementation and analytic gradients available in COLUMBUS allow for the characterization of relatively large molecules.

We have used COLUMBUS to characterize the excited states of *para*-benzyne,¹⁵⁵ and compared the geometries, energies, and electronic structures of the didehydro isomers of benzene,¹⁵⁶ pyrazine,¹⁵⁷ thiophene,¹⁵⁸ and furan. We have also utilized COLUMBUS to investigate aromatic tetraradicals formed through a double Bergman cyclization. The ability to perform geometry optimizations on molecules of this size has enabled us to report both adiabatic and vertical excitations, as well as analyze singlet–triplet gaps and through-bond and through-space couplings of the radical electrons. Additional features available in COLUMBUS utilized in these analyses included the determination of natural orbitals, unpaired electron densities, and spin densities.

We employed the state-averaging approach in COLUMBUS to demonstrate that *para*-benzyne exhibits a very dense manifold of excited valence and Rydberg states within a 10 eV window of the singlet 1Ag ground state.¹⁵⁵ In addition, spin-specific state averaging was instrumental in characterizing the 36 different spin states (1 quintet, 15 triplets, and 20 singlets) arising from 70 distinct electron configurations of a 1,4,5,8-tetra-dehydro-naphthalene. These results suggested a large (4 eV) gap between states 6 and 7 (Figure 17b), which
<figure><img src="image_27.png" /></figure>

<table><thead><tr><th></th><th>para<br>(ΔE(ST))</th><th>meta<br>(ΔE(ST))</th><th>ortho<br>(ΔE(ST))</th></tr></thead><tbody><tr><td>benzyne</td><td>49% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup> +<br>20% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup><br>5.99</td><td>62% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>7% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>21.76</td><td>65% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>4% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>38.21</td></tr><tr><td>pyrazine</td><td>58% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup> +<br>8% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup><br>28.27</td><td>58% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>11% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>-1.40</td><td>57% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>8% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>1.78</td></tr><tr><td>furan</td><td>59% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup> +<br>19% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup><br>7.83</td><td>60% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>13% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>19.35</td><td>58% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>15% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>15.96</td><td>63% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>10% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>14.54</td></tr><tr><td>thiophene</td><td>60% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>11% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>15.51</td><td>62% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>9% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>15.17</td><td>65% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>5% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>24.57</td><td>63% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> σ<sup>2</sup> +<br>8% π<sub>1</sub><sup>2</sup> π<sub>2</sub><sup>2</sup> π<sub>3</sub><sup>2</sup> (σ')<sup>2</sup><br>15.43</td></tr></tbody></table>

**Figure 17.** Examples of aromatic and heteroaromatic polyradicals characterized using COLUMBUS: (a) *para*-Benzyne contains a high density of valence and Rydberg excited states within 10 eV of the 1Ag ground state. Manifold determined using CAS (8,8) MRCI/TZ with 32-state averaging; (b) distribution of excited states for 1,4,5,8-tetradihydronaphthalene computing using state-specific, state-averaged CAS (6,6)/TZ geometry optimization; (c) dominant wave function configurations and singlet-triplet splittings (kcal/mol) for aromatic and heteroaromatic diradicals computed using singlet state MR-AQCC/TZ calculations and active spaces comprised of π and σ orbitals (CAS (8,8) for benzyne; CAS (8,7) for furan and thiophene, and CAS (12,10) for pyrazine including out-of-plane lone pair electrons on nitrogen atoms). Including the out-of-plane lone pair electrons on sulfur atoms in thiophene did not produce a significant difference, likely due to the mismatch in orbital energies.

then allowed us to focus our state-specific characterization efforts on the six lowest states.

We also utilized COLUMBUS to understand the multi-configurational character and singlet-triplet splitting in 9,10-didehydro-anthracene as well as in all didehydro isomers of pyrazine, benzyne, and thiophene (Figure 17c). We found that all these aromatic diradicals are two-configurational with wave functions comprising double occupation of the π orbitals and either the in-phase (σ) or out-of-phase (σ*) linear combination of the radical lobes. For the *ortho* and *meta* isomers, the dominant contribution to the wave function occupies the σ orbital. In contrast, for the *para* isomers, the σ*

component has the larger percentage for all systems except thiophene. The *para* isomers display the most significant through-bond coupling, while the *ortho* orientation typically results in the strongest through-space interaction. *meta* Isomers display the most distorted molecular geometries, likely to optimize through-space interactions in systems where through-bond coupling is minimal.

**4.5. Ionic Bonds in Authentic Hydrocarbons.** Since 1985, salts consisting solely of C and H have been isolated and characterized.¹⁵⁹⁻¹⁶² These results represent a breakthrough in the study of chemical bonds and indicate that steric crowding between the ions, along with an extended charge delocalization
<figure><img src="image_28.png" /><figcaption>**Figure 18.** CASSCF/MR-CISD results relevant for the formation of the ion pair between cyclopentadienyl anion and tropyl cation from the covalent compound: (a) Potential energy curves generated from single-point calculations (at the CASSCF level) on the geometries obtained from a relaxed scan along the C–C distance (shown in the structures) between the two fragments. The structure on the far right corresponds to the optimized geometry for the 3A' state. In both (b, c), the *f* and Δ*E*<sub>v</sub> results are given. In (b, c), they have been computed at the optimized geometries of the ground state and ion-pair, respectively. Both sets of Δ*E*<sub>v</sub> values (in eV) were calculated at the MR-CISD+Q level. The widths of the arrows represent relative oscillator strengths computed at the MR-CISD level.</figcaption></figure>

over at least one of the ions, are the two main factors that effectively suppress the formation of any covalent bond between the ions in the solid. In the limit of very extended charge delocalization over the two ions along with extreme steric congestion between the ions, one also has the suppression of any covalent bond in the solution.<sup>161,162</sup>

The study of ion pairs is crucial for understanding ionic compounds. In this context, some of the authors of this review recently studied a reduced model system,<sup>164</sup> which represents a good approximation to the steric repulsion between the tricyclopropylcyclopropenyl cation and Kuhn's anion studied in ref 160 It has been suggested through CASSCF and DFT calculations that a kinetically stable hydrocarbon ion pair can be formed from the corresponding radicals.<sup>164</sup> A kinetic model to study the time dependence of its concentration has also been developed.<sup>163</sup> DFT calculations indicate that the ion pair formed between the ions studied in ref 162 is significantly more stable than its covalent counterpart, a remarkable outcome for a hydrocarbon.<sup>165</sup>

Reactions between salts of the tropylium cation and the cyclopentadienyl anion usually form the covalent compound shown on the left of Figure 18a.<sup>166</sup> In this case, the parent ions are not sterically hindered enough and do not have extensive charge delocalization, which is consistent with the formation of the covalent compound. However, the potential energy curve for 3<sup>1</sup>A', shown in Figure 18a, suggests the formation of a hydrocarbon ion pair (containing two aromatic ions) in one of the excited states, as explained below.

The state-averaged CASSCF calculations performed with COLUMBUS provided a balanced description of several low-lying states important for the photodissociation of the C–C bond, as shown in Figure 18a. The potential energy curves shown in this plot were obtained from single-point calculations performed on the ground-state geometries generated from a relaxed scan (along the mentioned C–C bond) at the state-averaged CAS(12,11) level with the 6–31+G* basis set. C<sub>s</sub> symmetry was imposed, and five singlet states (3<sup>1</sup>A' + 2<sup>1</sup>A") were averaged with equal weights. At the optimized ground state geometry, the 11 active orbitals include a σ<sub>CC</sub>/σ<sub>CC</sub>* pair, five π and four π* orbitals, respectively. For larger distances, the σ<sub>CC</sub>/σ<sub>CC</sub>* pair turns into a π/π* pair, yielding six π and five π* orbitals, respectively.

Until ~2.5 Å, some states show large multiconfigurational characters involving complex electronic configurations. Further details will be published elsewhere. For larger distances, the configurations of four states (1<sup>1</sup>A', 2<sup>1</sup>A', 1<sup>1</sup>A", and 2<sup>1</sup>A"), due to the degeneracy of the SOMO orbitals of tropyl and cyclopentadienyl radicals, are consistent with dissociation yielding the mentioned ground state radicals. The fifth state (3<sup>1</sup>A') yields the tropylium cation and the cyclopentadienyl anion, both aromatic ions. The optimization (at the CASSCF level) of this state leads to the ion pair structure shown on the right side of Figure 18a. It is clear from the potential energy curves that, starting from the Franck–Condon (FC) region, this structure can be generated through nuclear relaxation either directly from the 3<sup>1</sup>A' state or via nonadiabatic transitions from other states.

Using the occupation numbers ($n_{occ}$) of the active orbitals at the CASSCF level, it was possible to set up an MR-CISD calculation yielding ~1.1 billion CSFs. The 6–31+G* basis set was adopted again. This calculation is based on a RAS + CAS + AUX scheme, where only single RAS → CAS and CAS → AUX excitations are allowed. The RAS and AUX orbitals were extracted from the initial set of 11 active orbitals, picking up the three with the largest and smallest values of $n_{occ}$, respectively. The remaining five orbitals were kept in the active space. This MR-CISD calculation was performed at the optimized geometries of the ground state (Figure 18b) and ion-pair (Figure 18c). Figure 18b shows the vertical excitation energies (Δ*E*<sub>v</sub>) at the MR-CISD+Q level, along with arrows whose widths represent the relative oscillator strengths (*f*).

The largest and the smallest *f* values are 0.200 and 0.004, corresponding to the 1<sup>1</sup>A' → 2<sup>1</sup>A" and 1<sup>1</sup>A' → 2<sup>1</sup>A' transitions, respectively (see Figure 18b). Thus, the former transition can be crucial in yielding the mentioned ion pair, depending on the potential energy curves at the MR-CISD/MR-CISD+Q level, which will be presented in a forthcoming publication. At this same level, the ion-pair energy is 3.47 eV larger than that of the ground state, as shown on the right side of Figure 18a. Assuming that the 1<sup>1</sup>A' → 2<sup>1</sup>A" transition can lead to the ion pair, we note that its Δ*E*<sub>v</sub> (6.32 eV, see Figure 18b) is slightly smaller than the approximate threshold energy of 6.49 eV required to generate the separated ions from the ground state structure. This value has been estimated by extrapolating the ground state potential energy curve to ~2.0
<figure><img src="image_29.png" /><figcaption>**Figure 19.** Selected active space orbitals (±0.04 isosurfaces) and occupations from CASSCF(6,20) included in the construction of the reference states of $U_2$ for the MRCI calculation of the $\Omega = 8_g$ state.</figcaption></figure>

eV (see Figure 18a) along with the adiabatic ionization energy of tropyl (6.28 eV)¹⁶⁷ and electron affinity of cyclopentadienyl (1.79 eV)¹⁶⁸ radicals. Similarly, the highest energy transition (¹¹A'→³¹A', see Figure 18b) is only 0.02 eV smaller than the mentioned threshold. Therefore, from an energetic point of view, all studied states in the FC region are capable of yielding the aromatic ion pair. However, as shown in Figure 18a, the profiles of the potential energy curves should also be taken into consideration. This point will be further discussed in a forthcoming publication.

Through comparison between Figure 18a and 18b it is clear that one has reordering of some states when dynamic electron correlation and extensivity correction are included. While at the CASSCF level the ordering is 2¹A' < 1¹A" < 3¹A' < 2¹A", at the MR-CISD+Q level it changes to 1¹A" < 2¹A' < 2¹A" < 3¹A'. It is also clear that the main effect is on the states of A" symmetry.

Assuming that the studied ion-pair is a minimum, a significant characteristic is its potential fluorescence. As shown in Figure 18c, it is a potentially fluorescent species in the near-infrared. The ion pair is also characterized by a significantly larger dipole moment of 15.47 D, compared to the dipole moment of 0.40 D in the ground state, both at the MR-CISD level. Its dipole moment is close to that of another hydrocarbon ion pair studied recently¹⁶⁴ and considerably larger than those of other ion pairs also studied at the MR-CISD level using COLUMBUS.¹⁶⁹⁻¹⁷³

**4.6. Uranium–Uranium Bond Analysis at the MRCI Level.** The COLUMBUS program package provides the basis for performing chemical bond analysis at the MRCI level of theory by utilizing the MRCI one-electron density matrix and the natural orbitals obtained from it. Transition metal and f-element dimers are known to be challenging systems for electronic structure methods, often requiring correlated multireference calculations. We exemplify the bond analysis here with the example of the uranium dimer, $U_2$, which has

previously attracted much attention and some controversy as to whether the bond order (BO) is 4 or 5.¹⁷⁴⁻¹⁷⁶

The calculations for $U_2$ were performed using the cc-pVTZ-DK3¹⁷⁷ basis set. CASSCF, restricted active-space SCF (RASSCF), and integral calculations were performed with OpenMolcas,²²,¹⁷⁸ using the second-order Douglas–Kroll–Hess Hamiltonian for treating scalar relativistic effects. Spin–orbit effects were not included. Bond analyses were performed with natural bond order (NBO) algorithms¹⁷⁹ using the NBO7 program¹⁸⁰ with input data generated with COLUMBUS. To render the NBO input generation convenient, a utility program, molcasto47, was recently developed by some of us.¹⁸¹

The interatomic distance of 2.56 Å was taken from the literature, based on the relativistic calculations for $U_2$ by Knecht et al.¹⁷⁶ A 6-electron 20-orbitals (6,20) active space was constructed following the work of Gagliardi and Roos.¹⁷⁴ The calculated occupations of the natural orbitals of the $\Omega = 8_g$ state (where $\Omega$ is the projection of the total angular momentum along the molecular axis) closely match the values reported in the literature,¹⁷⁴,¹⁷⁵ despite differences in the basis sets and bond parameters. Single-state CASSCF calculations identified the state with $\Omega = 8_g$ as the ground state, which was previously reported as one of the possible ground states.¹⁷⁴,¹⁷⁵,¹⁸² In subsequent RASSCF(12,26) calculations, additional $\sigma$ and $\pi$ bonding and antibonding orbitals were included with double excitations. Note that the ground state is studied here without consideration of the SO interaction and corresponds to the state $^7O_8$ with $\Lambda= 11$.

The MRCI calculations were based on the CASSCF-optimized orbitals. Here, 62 orbitals (including U shells up to 4s) were frozen, and a RAS(6,15) reference space was generated with just single excitations allowed to orbitals with CASSCF occupations smaller than 0.484. The active-space orbitals and allowed excitations were selected based on the occupations of CASSCF natural orbitals (Figure 19) of the $\Omega = 8_g$ ($B_{3g}$) state. The MRCI wave function included 22 configurations in the reference wave function following the
same state symmetry and allowed active space excitations. The selected reference state construction was based on an analysis of the CASSCF occupation numbers and set up with the intent to minimize the number of reference states and maximize the effects of dynamic correlation. All external double excitations conserving the state symmetry were allowed in the CI step for each reference state configuration; the total number of configurations in the MRCI calculation was approximately 100 million.

The effective bond order (EBO) definition used by Roos and Gagliardi¹⁷⁴ to assign U₂ a quintuple bond utilizes the typically noninteger occupation numbers obtained from multiconfigurational calculations. However, it does not explicitly consider the overlap of the participating atomic orbitals (AOs). The bonding and antibonding character of an orbital is qualitatively assessed to obtain the EBO, which is then defined as half of the sum of occupations of bonding valence MOs minus half of the sum of occupations of antibonding MOs. Based on the CASSCF occupations, the EBO is 4.14, similar to the value of 4.2 reported by Roos and Gagliardi. This result is assigned a quintuple bond character, considering that the bond order is above 4. Other numerical bond indices obtained in the present work from the 1-electron reduced density matrix and the AO overlap are collected in Table 3.

<table><caption>**Table 3.** U₂ Bond Orders Calculated with Different Wave Functions and Bond Order Criteria</caption><thead><tr><th>method</th><th>CASSCF(6,20)</th><th>RASSCF(12,26)</th><th>MRCI</th></tr></thead><tbody><tr><td>MBO</td><td>3.44</td><td>3.05</td><td>3.42</td></tr><tr><td>WBI</td><td>3.34</td><td>2.97</td><td>3.30</td></tr><tr><td>NLMO</td><td>2.01</td><td>4.01</td><td>4.01</td></tr></tbody></table>

The Wiberg bond index (WBI)¹⁸³ and Mayer bond order (MBO)¹⁸⁴ vary according to the level at which the dynamic correlation is treated and as a function of the active space. The WBI is defined for orthonormal orbitals, such as a set of Löwdin-orthonormalized atomic orbitals (AOs) or the natural atomic orbitals (NAOs) generated in NBO runs. Therefore, a measure of the AO overlap is implicitly included in the WBI. The MBO explicitly considers the AO overlap. We also report bond orders based on the natural localized molecular orbital (NLMO) population analysis generated in NBO runs; however, the results obtained for U₂ vary strongly and may therefore be unreliable. Note that bond indices that can be generated from the overlap and density matrices do not require a subjective assessment of whether a given MO is bonding or antibonding or whether it should be excluded from the BO estimation. They are, thus, arguably a more objective measure for bond order than EBO. As seen in Table 3, the WBI and MBO are similar for U₂ and considerably smaller than the EBO, indicating a BO of at most 4 for U₂. These findings agree with the BO assignment of U₂ by Knecht et al.¹⁷⁶

A critical factor in the assignment of the U₂ BO is the bonding furnished by the 5f AOs of $\phi$ symmetry. As shown in Figure 20, even with the relatively small isosurface values used for the plots, the resulting $\phi$-symmetric linear combinations are identified as barely bonding and only weakly antibonding, respectively. The MRCI orbitals appear similar to the $\phi$ orbitals of the active space. Their contributions to the covalent bonding in the dimer are, therefore, comparatively weak. Based on the occupations obtained from the different calculations

<figure><img src="image_30.png" /><figcaption>**Figure 20.** Bonding and antibonding $\phi$-symmetry natural localized molecular orbitals (from MRCI, ±0.01 isosurfaces).</figcaption></figure>

(Table 4), the contributions to the EBO are also small due to the nearly equal occupations of the bonding orbitals and their

<table><caption>**Table 4.** U₂ Occupations of $\phi$ and $\phi^*$ Symmetry Natural Localized Molecular Orbitals</caption><thead><tr><th>method</th><th>$\phi$</th><th>$\phi^*$</th></tr></thead><tbody><tr><td>CASSCF(6,20)</td><td>0.499</td><td>0.499</td></tr><tr><td>RASSCF(12,26)</td><td>0.496</td><td>0.494</td></tr><tr><td>MRCI</td><td>0.766</td><td>0.729</td></tr></tbody></table>

antibonding counterparts. However, the overall $\phi$ orbital populations are enhanced with MRCI, indicating that correlation outside the reference space influences the wave functions.

Overall, the present MRCI results indicate that the inclusion of dynamic correlation does not dramatically impact the resulting bond order of U₂ compared to what is obtained from multiconfigurational wave function calculations with carefully selected active spaces. However, the correlated MRCI results can be regarded as less active-space dependent and, therefore, somewhat more reliable because they include variationally calculated orbital populations across all orbitals beyond the frozen core–shells. More generally, therefore, MRCI is useful in this type of study, as it provides energies and wave functions at a consistently high level of theory, which is not possible for perturbative correlation methods. A more detailed analysis of the reference space structure and static correlations as they influence the numerical bond orders as well as the role of higher excitations is left for subsequent studies.

## 5. APPLICATIONS: EXCITATION ENERGIES AND POTENTIAL ENERGY SURFACES

### 5.1. Consistent Treatment of Covalent and Ionic States of Hexatriene at MR-CISD Level.
The study of all-trans polyene chains has a long-standing tradition due to their distinctive spectral properties, which arise from the structure of their low-lying excited-state manifold.⁵¹ In small- and medium-sized polyenes, the two low-lying excited states are the covalent $^1A_g^-$ and the ionic $^1B_u^+$ state. The wave functions of the covalent states are intricate, with significant contributions from doubly excited configurations, whereas those of the ionic states are predominantly composed of singly excited configurations.¹⁸⁵ This imposes significant challenges in accurately and uniformly describing both states.¹⁸⁵ MR-CISD and MR-AQCC calculations have previously been performed for ethylene¹⁸⁶ and butadiene,¹⁸⁷ the first members of this series. Here, we report an accurate and balanced treatment of both static and dynamic electron correlation effects in *trans* hexatriene using the variational MR-CISD and MR-AQCC methods, which are available in COLUMBUS. Our findings led us to extend our investigation to the electronic spectra of larger members of the all-*trans* polyene series.¹⁸⁸
A set of CASSCF molecular orbitals derived from a standard π-valence complete active space (CAS) with six electrons in six orbitals, CAS(6,6), was employed for the MR-CISD and MR-AQCC calculations. For these computations, the reference active space was constructed using the same π-valence CAS, and only the core orbitals were frozen. Vertical excitation energies from the ground state to the excited states $^{1}\text{B}_{\text{u}}^{+}$ and $^{1}\text{A}_{\text{g}^{-}}$ were calculated using the cc-pVnZ (n = D, T, Q) basis sets.⁶⁰ Extrapolation to the complete basis set (CBS) limit was performed using a two-point fit approach.¹⁸⁹ For the MR-CISD calculations, the Pople (+P) nonvariational size-extensivity correction was applied.¹⁸

Three main factors significantly influence the description of the two lowest excited states. These factors vary in magnitude depending on whether the state has covalent or ionic character. The first factor is the inclusion of dynamic σ-π electron correlation. At the CASSCF level, which primarily captures static correlation effects, the covalent state is well described, with no significant difference observed between the CASSCF and MR-CISD or MR-CISD+P results (Figure 21). In contrast,

<figure><img src="image_31.png" /><figcaption>**Figure 21.** CASSCF, MR-CISD, and MR-CISD+P vertical excitation energies (in eV) for the $^{1}\text{A}_{\text{g}^{-}}$ and $^{1}\text{B}_{\text{u}^{+}}$ states of hexatriene. The results were obtained using the cc-pVDZ basis set. The numerical difference between CASSCF and MR-CISD, as well as between MR-CISD and MR-CISD+P, are indicated for the $^{1}\text{B}_{\text{u}^{+}}$ state.</figcaption></figure>

the vertical excitation energy to the ionic state is strongly overestimated by 1.09 eV compared to the MR-CISD results. This overestimation increases to 1.81 eV when compared to the size-extensivity corrected MR-CISD+P result (Figure 21).

This discrepancy between the MR-CISD and MR-CISD+P results highlights the second factor influencing the description of the excited states: the size-extensivity error. While this error has a negligible impact on the covalent state, applying the size-extensivity correction stabilizes the ionic state by 0.72 eV (Figure 21). Notably, the MR-CISD+P results, corrected for size-extensivity, show good agreement with those obtained at the MR-AQCC level (Table 5).

The third factor is the basis set. For the covalent state, when comparing nonextrapolated excitation energies obtained using the cc-pVDZ basis set ($\Delta E_{\text{D}}$) with values extrapolated from double-$\zeta$ to triple-$\zeta$ ($\Delta E_{\text{DT}}^{\infty}$), a stabilization of 0.06 eV is

observed at the MR-CISD+P level (Table 5), increasing to 0.08 eV at the MR-AQCC level. For the ionic state, the stabilization amounts to 0.19 eV at the MR-CISD+P level and 0.27 eV at the MR-AQCC level (Table 5). Comparing $\Delta E_{\text{DT}}^{\infty}$ with excitation energies extrapolated from triple-$\zeta$ to quadruple-$\zeta$ ($\Delta E_{\text{TQ}}^{\infty}$), no significant difference is observed at the MR-CISD+P level for either state. However, at the MR-AQCC level, the ionic state undergoes further stabilization by 0.15 eV (Table 5). Notably, the extrapolated MR-AQCC results demonstrate slightly better numerical agreement with theoretical reference values.

Examination of the transition density (Figure 22) reveals that for the covalent state, the transition charges are centered along the bonds. In contrast, for the ionic state, they are localized on the atoms.¹⁷,³⁹ Moreover, the transition density of the ionic state at the MR-CISD level shows distinct σ-contributions as in-plane rings around the out-of-plane lobes that are absent in the CASSCF transition density. These σ-contributions are associated with the reduction of transition density self-repulsion, which is generally significant when there is considerable spatial overlap between the involved orbitals.⁴⁴ The studies by Kimber and Plasser⁴³,⁴⁴ on the role of transition density self-repulsion demonstrate that while the self-repulsion term raises the excitation energy, it can be compensated by $\sigma\sigma^*$ excitations. Consequently, the inclusion of σ-contributions in the wave function reduces the energetic penalty of the pure HOMO–LUMO transition, thereby lowering the excitation energy and resulting in reduced oscillator strengths.⁴,¹⁹²

In summary, while covalent states are well described by MCSCF theory, an accurate description of the ionic state requires extensive treatment of $\sigma-\pi$ electron correlation, the inclusion of size-extensivity corrections, and the use of a sufficiently large basis set. Accounting for these factors, the MR-CISD and MR-AQCC results consistently indicate that the ionic $^{1}\text{B}_{\text{u}^{+}}$ state is the first vertically excited state in trans-hexatriene. Based on the $\Delta E_{\text{TQ}}^{\infty}$ extrapolated results, this state lies 0.14 eV below the covalent state according to the MR-CISD+P calculations and 0.32 eV below according to the MR-AQCC calculations.

## 5.2. Charge-Transfer States in Doped Periacenes.
The generation of charge transfer (CT) states at donor/acceptor (D/A) interfaces is a crucial process that initiates hole/ electron separation in photovoltaic and photocatalytic processes.¹⁹³,¹⁹⁴ PAHs have gained significant interest as basic semiconductor materials due to their unique electronic band gap properties, which enable convenient tuning through doping. Doping with heteroatoms appears to be a promising way to provide flexible tuning possibilities. It has been shown previously¹⁹⁵⁻¹⁹⁹ that especially N or B,N doping is efficient in controlling the biradicaloid character of PAHs. In this work, borazine rings inserted into a (5a,5z) periacene sheet were used to function as a D/A interface. From several insertion

<table><caption>**Table 5.** MR-CISD+P and MR-AQCC Vertical Singlet Excitation Energies (in eV) of Hexatriene<sup>a</sup></caption><thead><tr><th rowspan="2">state</th><th colspan="3">MR-CISD+P</th><th colspan="3">MR-AQCC</th><th rowspan="2">TBE<sup>b</sup></th><th rowspan="2">CASPT2</th></tr><tr><th>ΔE<sub>D</sub></th><th>ΔE<sub>DT</sub><sup>∞</sup></th><th>ΔE<sub>TQ</sub><sup>∞</sup></th><th>ΔE<sub>D</sub></th><th>ΔE<sub>DT</sub><sup>∞</sup></th><th>ΔE<sub>TQ</sub><sup>∞</sup></th></tr></thead><tbody><tr><td><sup>1</sup>A<sub>g</sub><sup>-</sup></td><td>5.755</td><td>5.690</td><td>5.710</td><td>5.744</td><td>5.659</td><td>5.658</td><td>5.62</td><td>5.57<sup>c</sup></td></tr><tr><td><sup>1</sup>B<sub>u</sub><sup>+</sup></td><td>5.803</td><td>5.613</td><td>5.568</td><td>5.764</td><td>5.495</td><td>5.343</td><td>5.37</td><td>5.31<sup>d</sup></td></tr></tbody></table>

<footnote><sup>a</sup>Values computed using the cc-pVDZ basis set ($\Delta E_{\text{D}}$), as well as values extrapolated from double-zeta to triple-zeta ($\Delta E_{\text{DT}}^{\infty}$), and from triple-zeta to quadruple-zeta ($\Delta E_{\text{TQ}}^{\infty}$), are presented. Theoretical reference values are provided for comparison. <sup>b</sup>Theoretical best estimates (TBE) from ref 41. <sup>c</sup>Reference 190. <sup>d</sup>Reference 191.</footnote>
<figure><img src="image_32.png" /><figcaption>**Figure 22.** Transition densities (isovalue: 0.002 au) between the ground state and the excited states $^1A_g^-$ and $^1B_u^+$ of hexatriene, computed at the CASSCF and MR-CISD levels of theory using the cc-pVDZ basis set.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Figure 23.** (a) B,N-center and (b) B,N-vertical inserted borazine chains in a (5a,5z)-periacene. The boundaries of the borazine row are indicated by the black lines, which also define the segment numbers, marked by the red numbers.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 24.** Excitation energies (eV) and analysis of the character of the electronic transitions for (a) the BN-center and (b) BN-vertical structures. The state composition bar plot is numbered in terms of the three adopted segments (Figure 23), where the first number designates the hole and the second one the electron.</figcaption></figure>

possibilities investigated, the B,N-center, and B,N-vertical cases, as displayed in Figure 23, will be discussed in this contribution, as they exemplify well the creation of the CT states. Calculations on PAHs are challenging since they acquire open-shell character with an increasing number of benzene rings.<sup>11,200,201</sup> B,N doping is supposed to stabilize the PAH sheets, but the doped systems will possess a polyradical character as well.<sup>202</sup> Therefore, the use of multireference (MR)

calculations is crucial to obtain reliable results for ground and excited state properties.

The present contribution employs MR-CISD.<sup>50,51</sup> Energies include the *a posteriori* Pople size extensivity correction,<sup>18</sup> which is denoted as MR-CISD+P. Molecular orbitals used in the MR-CISD calculation are derived from a state-averaged (SA)-MCSCF calculation. The reference space of the MR-CISD and the active space of the MCSCF calculation use the complete active space (CAS) of eight electrons in eight orbitals
CAS(8,8). The eight π orbitals, corresponding to the b₁ and a₂ irreducible representations, were chosen according to natural orbital (NO) occupations. In each of the two irreducible representations, the two highest occupied (HONOs) and the two lowest unoccupied (weakly occupied) NOs (LUNOs) were chosen. The choice of this reference space has been made based on the weight of the reference configurations in the wave function expansion (nonreference configurations should have a weight of less than approximately 2%), on NO occupations, and experience obtained with previous calculations on the typical biradical compound heptazethrene.¹⁴⁷ For reasons of computational efficiency, all σ orbitals were frozen. It has also been previously shown that for the calculation of energy differences between singlet/triplet or different spin multiplicities for phenalenyl, freezing all σ orbitals influences the results by ~0.1 eV only.¹⁰⁰,¹⁴⁷ Moreover, the 6-31G basis set¹⁴⁸ was used following the work reported in ref 147 where differences of <0.1 eV were found for excited states of phenalenyl in comparison to larger, polarized basis sets. The wave function analysis was carried out through THEODORE using the one-particle transition density matrices.⁴⁰,²⁰³ More details can be found in ref 204.

Figure 24 presents the energy diagram of the four lowest excited states for the BN-center and BN-vertical doped periacene structures, classifying them as locally excited (LE) or charge-transfer (CT) states. The excitations are decomposed among the three molecular fragments defined in Figure 23. In the case of the BN-center structure, the bright state is S₂, which is dominated by a “3-3” type state, characterized by a locally excited configuration on the bottom pentacene. S₁ exhibits a significant “3-1” type CT character, meaning an excitation from the bottom pentacene to the top pentacene unit. S₃ is a state completely dominated by double excitations (as defined in ref 205). It is, therefore, not accessible by single-excitation methods such as the popular TDDFT method or the second-order algebraic construction method (ADC(2)). It is located above the bright state but is quite close in energy, so one might want to consider it as well. The BN-vertical structure describes an even more critical situation concerning the doubly excited state: it is S₁ lying well below 1 eV and, therefore, all other states. The bright state is S₃. Reviewing Figure 24, we note that aside from the S₃ and S₁ excited states discussed above, the different states also have a significant admixture of doubly excited character (indicated by light yellow bars representing the Ω-diagnostic²⁰⁵), highlighting the challenges of accessing any of these states without a multireference framework.

Natural transition orbitals (NTOs) for the states dominated by single excitations and electron density differences for the doubly excited states are shown in Figure 25 for selected transitions. The S₁ state of the BN-center structure represents a conventional CT state, which is separated by the borazine interface (as discussed above). It should be well-suited for subsequent charge separation dynamics. Interestingly, the two mentioned doubly excited states exhibit a similar CT character, where the borazine chain also serves as a boundary.

Both BN-center and BN-vertical structures present interesting CT cases and potentially could be used as photocatalysts in oxidation/reduction processes, given the separation of the electron/hole pair. We note that, especially for the BN-vertical doped structure, consideration of double excitations is essential. Therefore, this work demonstrates the need for

<figure><img src="image_35.png" /><figcaption>**Figure 25.** NTOs showing hole and electron orbitals and density difference plots between the excited and ground states for (a) the BN-center and (b) the BN-vertical structure. For the density difference plots, the red color indicates a decrease in electron density, and blue indicates an increase in density. H stands for highest occupied MO (HOMO) and L for lowest unoccupied MO (LUMO).</figcaption></figure>

high-level computational methods to accurately understand the electronic states of these doped periacene structures.

### 5.3. Spin-Orbit CI on the Electronic Structure of the CUO Molecule.
The Spin-Orbit Configuration Interaction (SO-CI) formulated with GUGA¹³ represents a robust computational methodology in quantum chemistry, particularly for systems where both electron correlation and spin-orbit coupling play a significant role. In the COLUMBUS GUGA SO-CI, the spin-orbit coupling terms are incorporated directly into the CI Hamiltonian, allowing for the simultaneous treatment of electron correlation and relativistic spin-orbit effects. The COLUMBUS programs initially implemented the relativistic effective core potentials²⁰⁶ and the associated spin-orbit operators.²⁰⁷ Subsequently, the implementation was adapted to work with all-electron relativistic Hamiltonians, using mean-field spin-orbit integrals as available in OpenMolCAS.¹⁴,²⁴,²⁰⁸ Furthermore, the spin-orbit operators and integrals implemented in the COLUMBUS programs, as well as the spin-orbit configuration interaction, can also be used in the implementation of quantum computing for Hamiltonians with spin-orbit interactions, for example, refs 209,210.

Spin-orbit coupling is crucial for describing heavy-element systems, where these effects have a significant impact on spectroscopic and chemical properties. The GUGA CI with spin-orbit interaction provides a robust framework for understanding complex phenomena, such as fine-structure splitting, spin-forbidden transitions, and the electronic structure of systems with significant relativistic contributions. A wide range of chemically and spectroscopically relevant systems has been studied.²¹¹⁻²¹⁵

In this work, we review the application of the GUGA CI with spin-orbit interactions to the CUO molecule and illustrate the approach's power for accurate calculations, as well as its capability to resolve discrepancies in the extensive experimental and theoretical work. The spectroscopic properties of the molecule CUO have been extensively studied to establish the structural properties and low-lying electronic states of the CUO molecule. A study of the infrared spectra of CUO in rare-gas matrices revealed that the molecule is linear and that the two stretching frequencies differed by unusual
amounts when the host gas was changed from neon (Ne) to argon (Ar).<sup>216,217</sup> The study established the interaction of inert gases with the polar molecules containing uranium atoms. Extensive calculations have been carried out to understand the electronic structures of the ground and low-lying electronic states, as well as the weak interaction with the noble gas matrix, employing a range of methods from DFT calculations to sophisticated coupled cluster calculations while treating spin-orbit effects in various approximations (see ref 215 and references therein). However, it is only with the spin-orbit interaction GUGA CI approach that the correct order of the low-lying electronic states is established, indicating that large-scale configuration interactions with correlation and spin-orbit effects treated on the same footing are necessary for an accurate determination of the electronic structures in such cases.

Figure 26 shows the low-lying orbitals and electronic states, as well as the CU and OU stretch frequencies in the low-lying

<table><caption>**Figure 26.** Electronic structure of CUO. The CU and OU vibrational frequencies are the experimental measurements of electronic ground states. In the Ne matrix, when the interaction with the inert gas is minimal, the closed-shell state is the ground state. In the Ar matrix, the open-shell polarized electronic state becomes the lower-lying electronic state due to the strong interaction of Ar with the polarized state. The pictures shown here depict the electronic configurations at the SCF level of calculations.</caption><thead><tr><th></th><th></th><th></th><th>Low lying states:</th></tr></thead><tbody><tr><td>— —</td><td>1φ U f</td><td></td><td>8σ², <sup>1</sup>Σ, Ω=0</td></tr><tr><td>— —</td><td>1δ U f</td><td></td><td>ν<sub>CU</sub> = 1047 cm<sup>-1</sup></td></tr><tr><td>— —</td><td></td><td></td><td>ν<sub>OU</sub> = 872 cm<sup>-1</sup></td></tr><tr><td>— —</td><td>4π C p + U f,d</td><td></td><td>8σ<sup>1</sup>φ<sup>1</sup>, <sup>3</sup>φ, Ω=2,3,4</td></tr><tr><td>— —</td><td>8σ C p + U f,d</td><td></td><td>ν<sub>CU</sub> = 853 cm<sup>-1</sup></td></tr><tr><td>— —</td><td></td><td></td><td>ν<sub>OU</sub> = 804 cm<sup>-1</sup></td></tr><tr><td>— —</td><td>7σ O p + U f,d</td><td></td><td>8σ<sup>1</sup>1δ<sup>1</sup>, <sup>3</sup>Δ, Ω=1,2,3</td></tr><tr><td>— —</td><td>3π O p + U f,d</td><td></td><td></td></tr></tbody></table>

electronic states. Experimentally, in Ne matrices, the infrared spectra indicate that the closed shell state with the $\sigma^2$ electronic configuration is the ground state. However, in the Ar matrix, the order of the low-lying electronic states changed due to the larger interaction of Ar with the open shell triplet states with electronic configuration $\sigma^1\varphi^1$. DFT, CASMP2, and CCSD(T) calculations all place the triplet state as the ground state in the absence of spin-orbit interactions and position the triplet state further below the closed-shell state when spin-orbit interactions are considered. The effect of the Ar interaction is estimated to be 0.05 eV; thus, those calculations failed to explain the observed reversal of the Ar interaction's effect on the order of low-lying electronic structures compared with that in the Ne matrix. Our extensive spin-orbit CI calculations show that when spin-orbit interaction and correlation effects are considered simultaneously, the ground state of CUO is calculated to be the singlet state of the closed-shell configuration, $^1\Sigma$, $\Omega = 0$. The lowest electronic state from the open-shell triplet configuration is calculated to be the $\Omega = 2$ state, which is predominantly $^3\Delta$ of the $\sigma^1\delta^1$ configuration. It is only ca. 0.047 eV above the calculated ground state energy. Our results are consistent with the experimental observation that the closed-shell singlet state is the ground state observed in the Ne matrix. The stronger interaction with the Ar matrix, estimated to be 0.005 eV, can

lower the triplet configuration state to become the ground state in the Ar matrix.

# 6. APPLICATIONS: PHOTODYNAMICS

## 6.1. Nonadiabatic Molecular Dynamics with the COLUMBUS/Newton-X Interface.
COLUMBUS has often been used in nonadiabatic molecular dynamics simulations.<sup>218-220</sup> Its analytical gradients and nonadiabatic couplings are invaluable for this type of application. Indeed, when the trajectory surface hopping platform Newton-X was developed, its first interface was with COLUMBUS.<sup>217</sup> It provided information on electronic structures, such as potential energies, energy gradients, and nonadiabatic couplings for dynamics. COLUMBUS still plays a crucial role in Newton-X, which has recently extended its capabilities to run ab initio multiple spawning (AIMS) using the methods available in COLUMBUS.<sup>218</sup>

Trajectory surface hopping propagates a swarm of independent classical trajectories.<sup>221</sup> At each integration time step, a trajectory lies on a single Born-Oppenheimer potential energy surface of one of the molecular electronic states. To emulate the nonadiabaticity of the nuclear wavepacket (i.e., the population transfer between states), surface hopping estimates transition probabilities between pairs of states and stochastically allows the classical trajectories to switch to another state. Thus, the ensemble of independent trajectories forms a nonadiabatic classical wavepacket by capturing branching events.

The propagation of independent trajectories and the need for electronic structure information exclusively at classical nuclear geometries make surface hopping suitable for on-the-fly (or direct) propagation,<sup>222</sup> where the electronic structure is computed during the dynamics instead of being stored as a precomputed multidimensional potential energy surface.

Nonadiabatic dynamics with COLUMBUS can be directly run using its native nonadiabatic coupling vectors, time-derivative couplings based on wave function overlaps,<sup>223,224</sup> or still curvature-based time-derivative couplings using the time-dependent Baeck-An approach.<sup>225</sup> It can also be done on complex-valued potential energy surfaces.<sup>226</sup>

One feature distinguishing COLUMBUS from most other multiconfigurational quantum chemistry programs is its outstanding flexibility in building customized wave functions. Thus, when faced with the challenge of describing excitations and processes that would require too large active orbital spaces, such as dissociation or bond formation, COLUMBUS can provide ways to reduce this description by any creative arrangement of subspaces and electron transitions between them.

This flexibility is illustrated in the surface hopping dynamics of ethylene reported in ref 219 in that study, the ultrafast dynamics of the $\pi\pi^*-\text{excited}$ ethylene were simulated, including Rydberg states and the possibility of dissociation. Accounting for those effects required a complete active space of 12 electrons in 16 orbitals, which is too computationally expensive for on-the-fly evaluations during surface hopping. Nevertheless, with COLUMBUS, the following combination of subspaces was built at the MCSCF level
<figure><img src="image_36.png" /><figcaption>Figure 27. Molecular orbitals in the generalized active space used in the surface hopping of cyclobutanone. Reproduced with permission from ref 220. Copyright 2024 American Institute of Physics.</figcaption></figure>

$$ \begin{aligned} \text{MCSCF} = & \{[\text{PP}_{\text{CC}}] \times 4[\text{PP}_{\text{CH}}] \times [\text{CAS}(2, 2)_{\pi} \\ & \rightarrow 4\text{AUX}_{\text{Ryd}}]\} \\ & + \{[\text{PP}_{\text{CC}}] \times 3[\text{PP}_{\text{CH}}] \times [\text{CAS}(2, 2)_{\pi} \leftrightarrow \text{PP}_{\text{CH1}}]\} \\ & + \{[\text{PP}_{\text{CC}}] \times 3[\text{PP}_{\text{CH}}] \times [\text{CAS}(2, 2)_{\pi} \leftrightarrow \text{PP}_{\text{CH2}}]\} \\ & + \{[\text{PP}_{\text{CC}}] \times 3[\text{PP}_{\text{CH}}] \times [\text{CAS}(2, 2)_{\pi} \leftrightarrow \text{PP}_{\text{CH3}}]\} \\ & + \{[\text{PP}_{\text{CC}}] \times 3[\text{PP}_{\text{CH}}] \times [\text{CAS}(2, 2)_{\pi} \leftrightarrow \text{PP}_{\text{CH4}}]\} \end{aligned} \quad (34) $$

In this partition, $PP_{AB}$ consists of a subspace containing two orbitals ($\sigma_{AB}$ and $\sigma_{AB}^*$) and two electrons each, where only the two doubly occupied configurations are allowed. The wave function also allowed single excitations from the $CAS(2,2, \pi, \pi^*)$ space into four Rydberg orbitals. Finally, single and double substitutions from the $CAS(2,2)$ into any $PP_{AB}$ space (and vice versa) were included.

Using this MCSCF wave function, MR-CISD was computed, including CSFs arising from single and double excitations from the following reference space

$$ \begin{aligned} \text{MRCI}_{\text{ref}} = & [\text{CAS}(2, 2)_{\pi} \\ & \rightarrow 4\text{AUX}_{\text{Ryd}}] + [\text{CAS}(4, 4)_{\pi\sigma(\text{CH1})} \\ & + [\text{CAS}(4, 4)_{\pi\sigma(\text{CH2})}] + [\text{CAS}(4, 4)_{\pi\sigma(\text{CH3})}] \\ & + [\text{CAS}(4, 4)_{\pi\sigma(\text{CH4})}] \end{aligned} \quad (35) $$

More recently, COLUMBUS' wave function flexibility was explored to describe the nonadiabatic dynamics of cyclobutanone excited into a Rydberg state, which could induce diverse types of dissociation.²²⁰ In this case, the MCSCF wave function was chosen as

$$ \text{MCSCF} = 4[\text{PP}_{\text{CC}}] \times [\text{PP}_{\text{CO}}] \times [\text{CAS}(4, 4)_{n, \pi, \pi^*, 3s}] \quad (36) $$

with the molecular orbitals shown in Figure 27.

The construction of such specialized wave functions requires a high level of conceptual and technical expertise. Nevertheless, it pays off in terms of the superior quality of the results, with the flexibility to describe the photochemical process at

relatively low costs, which cannot be achieved with standard active-space procedures.

Beyond the direct use of COLUMBUS in nonadiabatic dynamics, the advent of machine learning (ML) techniques for computational chemistry has opened other potential usages for the program. In this case, COLUMBUS is employed as a quantum chemistry program to build high-level data sets for training machine learning potentials. Such a novel application is illustrated in ref 227, where COLUMBUS was used during active learning of an ML potential for fulvene. The converged ML potential required about 6000 single-point evaluations of energies and energy gradients at the CASSCF level.

The advantage of such a procedure is that it significantly reduces computational costs. In the fulvene case, the computational cost due to the 6000 COLUMBUS calculations is equivalent to running only ten surface hopping trajectories. Nevertheless, after converging the ML potential, we can run thousands of new trajectories at a minimal cost.

## 6.2. Nonadiabatic Dynamics Using Machine Learning in SHARC.
One of the first ML studies aimed at accelerating nonadiabatic surface hopping molecular dynamics utilized a reference data set derived from MR-CISD calculations performed with the COLUMBUS software package.²²⁸ The exceptional accuracy of electronically excited states offered by MR-CISD, combined with the ability to calculate nonadiabatic coupling vectors, made COLUMBUS an ideal tool for this application.

The goal here was to facilitate the interpolation of quantum chemical properties at significantly lower computational costs while maintaining the accuracy of advanced quantum chemical methods. To achieve this, an analytical function must be created that maps input variables, usually molecular geometry (X), to output properties (Y), such as potential energy surfaces, gradients, dipole moment vectors, and nonadiabatic coupling vectors. Unlike traditional physical models, which rely on explicit assumptions about the system, machine learning models utilize universal approximators that can represent any function with arbitrary accuracy, provided sufficient training data and model complexity are available. Consequently, ML proves to be highly effective when integrated into surface
<figure><img src="image_37.png" /><figcaption>**Figure 28.** Concept of ML for surface hopping molecular dynamics, including training of ML potentials and training set generation. Adapted with permission from ref 228. Available under a CC-BY license. Copyright Westermayr et al.</figcaption></figure>

hopping molecular dynamics, provided there is suitable training data and mechanisms to monitor the reliability of the ML-generated potentials during simulations.<sup>229–231</sup> These principles have been integrated into the SHARC program (surface hopping including arbitrary couplings),<sup>232–234</sup> which pioneered the propagation of trajectories concurrently driven by nonadiabatic couplings and spin-orbit couplings.<sup>235</sup>

The SHARC method interfaces with COLUMBUS, as well as with many other quantum chemistry programs, which can provide energies, gradients, and ideally nonadiabatic couplings and spin-orbit couplings. In the absence of the latter, SHARC will evaluate these numerically, thereby uniquely positioning itself to handle internal conversion and intersystem crossing on an equal footing. Moreover, SHARC can be integrated with neural networks, such as SchNarc,<sup>236</sup> which is based on the deep continuous filter of the convolutional neural network SchNet.<sup>237,238</sup> It can also be integrated with the recently updated SPaiNN,<sup>239</sup> which uses the equivariant polarizable neural network for atomic interactions PaiNN,<sup>240</sup> an equivariant version of SchNet. This integrated framework, as illustrated in Figure 28, enables the simulation of molecular dynamics on a long time scale based on training data from COLUMBUS.

ML modules can utilize the SHARC tools, such as Wigner sampling, to generate initial training sets for ML models, which are then employed to train neural network potentials for surface hopping. Additionally, a technique known as active learning or adaptive sampling is used to assess the reliability of the machine learning predictions. This approach ensures continuous refinement of the ML models in areas where their predictions deviate from quantum chemical reference calculations. This is achieved by comparing multiple ML

predictions for each new molecular geometry accessed during photodynamics. If discrepancies arise between the ML models, quantum chemical reference calculations are performed, and the ML models are then retrained to restore accuracy. This iterative process enables longer simulation times as the accuracy of the ML potentials is maintained throughout the dynamics. These tools have been leveraged in SHARC to perform nanosecond-scale simulations for the methylenimmonium cation, the smallest member of protonated Schiff bases used to model the processes involved in vision.<sup>228</sup> As the main chemical reaction is a cis-trans isomerization, multireference methods, such as those implemented in COLUMBUS, are essential for accuracy; however, propagating trajectories for over a nanosecond using COLUMBUS on the fly alone would not have been feasible.

Since this initial exploration of long-time scale photodynamics using ML, several additional methods and studies have been developed, including MLAtom 2 and PyR-AI²MD.<sup>241,242</sup> Recently, photodynamics simulations of the amino acid tyrosine were performed in conjunction with SHARC by integrating training data from various quantum chemical reference methods. This integration facilitated the discovery of roaming atoms—radicals that migrate around molecules and potentially lead to photodamage.<sup>243</sup> While these studies have revealed intriguing aspects of natural photodynamic processes, numerous challenges remain in extending ML-driven photodynamics to more complex systems, particularly those with environmental influences.

One of the biggest challenges is the transferability of ML potentials across chemical compound space.<sup>244,245</sup> Generating reference data for excited states is much more difficult.<sup>229</sup> and accurate modeling of conical intersections—regions where
<figure><img src="image_38.png" /><figcaption>**Figure 29.** S₀, T₁, and S₁ energy levels (relative to S₀ of structure A) at MR-CISD and MR-AQCC (upper panel, MRAQCC energies in parentheses), and Mukherjee's MR-CCSD(T) levels (bottom panel) with CAS(2,2) active space in the aug-cc-pVDZ basis set for several relevant structures. The geometries (except for E) have been optimized at the SA-CASSCF(2,2)/aug-cc-pVDZ level as follows: (A) Optimized S₀, (B) optimized T₁, (C) optimized S₀ with constraints, (D) optimized T₁ with constraints, (E) geometry from MD, (F) optimized S₀ of SiH₃CH diradical, (G) optimized T₁ of SiH₃CH diradical, (H) optimized S₀ of CH₃SiH diradical, (I) optimized T₁ of CH₃SiH diradical. All structures except C–E are local minima.</figcaption></figure>

electronic states converge and facilitate nonradiative transitions—is essential but challenging to achieve with ML models, which are typically smooth functions. Recent advances in the use of companion matrices show promise in solving this problem.²⁴⁴,²⁴⁵ Despite these obstacles, ML-driven dynamics, as demonstrated in studies of molecules such as the methylenimmonium cation, show that neural networks can provide both efficiency and precision for complex processes in excited states. Furthermore, high-quality data obtained with programs such as COLUMBUS are particularly beneficial for these fundamental studies.

**6.3. Toward MD of Silaethylene with Nonadiabatic and Spin-Orbit Couplings.** Molecules with photoexcited π-bonds exhibit fascinating dynamics that involve a complex interplay of torsion, pyramidalization, and hydrogen transfer en route to the conical intersection seam, facilitating a transition to the electronic ground state. Silaethylene (SiCH₄) is a prime example of a polar π-bond and a computationally feasible model for ab initio nonadiabatic molecular dynamics (NAMD)

studies. The previous NAMD study of silaethylene by Zechmann et al.²⁴⁶ focused on low-lying singlet potential energy surfaces. However, incorporating intersystem crossing (ISC) via the inclusion of the triplet manifold into the analysis is crucial, as ISC can be significant in systems with non-negligible spin-orbit coupling due to the presence of heavier atoms.

The availability of analytical gradients, nonadiabatic couplings, and spin-orbit couplings at the MRCI level within the COLUMBUS program makes this code particularly suitable for accurate NAMD studies of smaller molecules.¹³,¹⁴,²⁴⁷,²⁴⁸ In previous studies, we implemented surface hopping MD with nonadiabatic and spin-orbit couplings by interfacing the Newton-X program with different quantum chemistry codes, and presently, we extended this interface to COLUMBUS, including spin-orbit couplings from MRCI.²⁴⁹⁻²⁵¹

However, it is crucial to address a well-known limitation of the MRCI method: the lack of size extensivity, which can lead
to inaccuracies in energy calculations. To verify the reliability of the MR-CISD results, we compare them to the MR-AQCC and Mukherjee's multireference coupled cluster theory, incorporating triples perturbatively (MR-CCSD(T)).<sup>25,253</sup> While MR-AQCC gradients have been implemented,<sup>25,5</sup> there are no NAD and SO couplings available for either method, meaning that they cannot be used in the MD directly.

We employed the SA-CASSCF (2,2) method in the aug-cc-pVDZ basis to optimize different structures of silaethylene in their S₀ and T₁ states. We have investigated the planar silaethylene structure, as well as structures with a diradical center at carbon and silicon. We also investigated a 90-degree twisted structure and one structure along an MD trajectory, where the T₁ state was below the S₀. For these geometries, we then performed single-point calculations at the MR-CISD, MR-AQCC, and Mukhejee's MR-CCSD(T) levels. We computed three lowest-lying states (S₀, S₁, and T₁) for MR-CISD and two states (S₀ and T₁) for MR-AQCC using COLUMBUS. Mukherjee's MRCCSD(T) calculations have been performed using Orca.<sup>254</sup>

The results of the aforementioned methods for the nine structures of SiCH₄ (labeled A to I) are presented in Figure 29. Structures A and B represent the optimized S₀ and T₁ local minima of silaethylene, respectively. Structure A represents the global minimum, serving as the reference for the relative energies displayed. Geometries C and D were optimized with the Si–C dihedral angle constrained to 90° and are not local minima. A twisted structure with T₁ lying below S₀ at the MR-CISD level was obtained from a molecular dynamics trajectory and is denoted by E. Finally, we considered two diradical isomers of silaethylene: SiH₃CH and CH₃SiH. F and G denote the S₀ and T₁ optimized structures of SiH₃CH, while H and I represent the S₀ and T₁ optimized structures of CH₃SiH. Geometries F–I have been verified as local minima in the respective electronic state. For MR-CISD, we also computed the S₁ energies.

The relative energies of all structures, as well as the separations between the singlet and triplet states, show an excellent agreement across all three computational methods. For structure B, we also computed the MkCCSD(T) S-T gap in the aug-cc-pVTZ basis, yielding 1.1 kcal/mol; in contrast, the results in the aug-cc-pVDZ basis are 0.65 kcal/mol. The aug-cc-pVDZ basis, which is computationally affordable for the MD study, thus appears to be a suitable choice.

In some cases, when the energy separation between singlet and triplet states is minimal, the relative order of the two states changes depending on the method. This happens particularly for the nonequilibrium structures C, D, and E, which are easily energetically accessible in an MD started from S₁ at the equilibrium geometry A. In our opinion, the swapping of extremely closely lying states should not have an adverse influence on the MD. The diradical equilibrium structures F–I have a relatively weak multireference character (second largest CI coefficient about 0.2) and have larger singlet–triplet gaps than the nonequilibrium structures C-E. The SiH₃CH diradical (F, G) always has a triplet ground state and is relatively high in energy compared to the stable molecule A. However, the CH₃SiH diradical structure (H and I) has a singlet ground state. Its singlet local minimum (H) is almost degenerate in energy to the silaethylene molecule A and, depending on the energy barrier might be a stable isomer and might be reachable as a product of a photoisomerization of the molecule A after excitation to S₁, which might also include triplet intermediates.

A photorelaxation back to A is another option, as is the generation of products in the triplet state, which is energetically accessible. A NAMD study, including the SOC effects, will be needed to elucidate these processes.

To summarize, our results confirm that MR-CISD is both computationally efficient and reliable as a basis for NAMD simulations of the photochemistry of silaethylene after excitation to the S₁ state, which will be the subject of a forthcoming paper.

## 7. CONCLUSIONS AND OUTLOOK

The topics covered in this work illustrate the broad applicability of COLUMBUS in recent years. The need for multireference treatments is particularly evident in the section on open-shell systems, where di- and polyradicals serve as key examples. In the context of excited states, we addressed the challenge of achieving a balanced description of ionic and covalent states in hexatriene, as well as the construction of charge-transfer and doubly excited states. Spin-orbit MRCI calculations exemplify the extension of COLUMBUS into more demanding theoretical areas.

COLUMBUS has long been at the forefront of surface-hopping dynamics, owing to the availability of analytic energy gradients and nonadiabatic coupling vectors. However, the computational cost of performing MRCI-based on-the-fly dynamics has traditionally limited its application to small systems. As shown in the photodynamics section, this limitation is now being overcome by integrating machine learning techniques, which promise to significantly extend the scope of these simulations.

Recent code developments span several areas, including the computation of spin densities at the MCSCF level and the implementation of empirical scaling procedures designed to correct the imbalance between ionic and covalent states. While this issue has also been addressed by systematically expanding the wave function, scaling approaches provide a more computationally efficient alternative. The treatment of electronic resonances using complex absorbing potentials and simulations of the dynamics of highly excited states involved in attosecond photoemission spectroscopy further demonstrate COLUMBUS's ability to tackle problems beyond the scope of conventional quantum chemistry.

In summary, the diverse topics presented here are unified by a common foundation in multireference theory. Together, they highlight both the feasibility and the necessity of applying multireference methods to an increasingly diverse range of complex chemical and spectroscopic problems.

### AUTHOR INFORMATION

#### Corresponding Authors

Felix Plasser – Department of Chemistry, Loughborough University, Loughborough LE11 3TU, United Kingdom; ● orcid.org/0000-0003-0751-148X; Email: f.plasser@lboro.ac.uk

Hans Lischka – Department of Chemistry and Biochemistry, Texas Tech University, Lubbock, Texas 79409, United States; ● orcid.org/0000-0002-5656-3975; Email: Hans.Lischka@ttu.edu

#### Authors

Ron Shepard – Chemical Sciences and Engineering Division, Argonne National Laboratory, Lemont, Illinois 60439, United States; ● orcid.org/0000-0002-0272-0824
Péter G. Szalay — Laboratoy of Theoretical Chemistry,
Institute of Chemistry, ELTE Eötvös Loránd University, H-
1117 Budapest, Hungary; ● orcid.org/0000-0003-1885-
3557

Russell M. Pitzer — Department of Chemistry and
Biochemistry, The Ohio State University, Columbus, Ohio
43210, United States

Rodolpho L. R. Alves — Departamento de Química, CCEN,
Universidade Federal da Paraíba, 58059-900 João Pessoa,
Brazil

Adelia J. A. Aquino — Department of Mechanical Engineering,
Texas Tech University, Lubbock, Texas 79409, United
States; ● orcid.org/0000-0003-4891-6512

Jochen Autschbach — Department of Chemistry, University at
Buffalo, State University of New York, Buffalo, New York
14260-3000, United States; ● orcid.org/0000-0001-9392-
877X

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

Jhonatas R. Carvalho — Department of Chemistry and
Biochemistry, Texas Tech University, Lubbock, Texas 79409,
United States

Julio C. V. Chagas — Department of Chemistry and Advanced
Scientific Computing and Modeling Laboratory, Aeronautics
Institute of Technology, 12228-900 São José dos Campos,
Brazil; Department of Chemistry, Northwestern University,
Evanston, Illinois 60208, United States; ● orcid.org/0000-
0001-9647-0821

Leticia González — Institute of Theoretical Chemistry, Faculty
of Chemistry and Vienna Research Platform on Accelerating
Photoreaction Discovery, University of Vienna, 1090 Vienna,
Austria; ● orcid.org/0000-0001-5112-794X

Andreas Hansen — Mulliken Center for Theoretical Chemistry,
University of Bonn, 53115 Bonn, Germany; ● orcid.org/
0000-0003-1659-8206

Bhumika Jayee — Department of Chemistry and Biochemistry,
Texas Tech University, Lubbock, Texas 79409, United
States; Physical Sciences and Computational Division, Pacific
Northwest National Laboratory, Richland, Washington
99354, United States

Miklos Kertesz — Chemistry Department and Institute of Soft
Matter, Georgetown University, Washington, District of
Columbia 20057-1227, United States; ● orcid.org/0000-
0002-7930-3260

Francisco B. C. Machado — Department of Chemistry and
Advanced Scientific Computing and Modeling Laboratory,
Aeronautics Institute of Technology, 12228-900 São José dos
Campos, Brazil; ● orcid.org/0000-0002-2064-3463

Spiridoula Matsika — Department of Chemistry, Temple
University, Philadelphia, Pennsylvania 19122, United States;
● orcid.org/0000-0003-2773-3979

Silmar A. do Monte — Departamento de Química, CCEN,
Universidade Federal da Paraíba, 58059-900 João Pessoa,
Brazil

Saikat Mukherjee — Faculty of Chemistry, Nicolaus
Copernicus University in Torun, 87100 Torun, Poland;
● orcid.org/0000-0002-0025-4735

Dana Nachtigallová — Institute of Organic Chemistry and
Biochemistry, Academy of Sciences of the Czech Republic, CZ-
16610 Prague, Czech Republic; IT4Innovations,
VŠBTechnical University of Ostrava, 708 00 Ostrava, Czech
Republic; ● orcid.org/0000-0002-9588-8625

Reed Nieman — Department of Chemistry and Biochemistry,
Texas Tech University, Lubbock, Texas 79409, United
States; ● orcid.org/0000-0003-3011-5223

Vytor P. Oliveira — Department of Chemistry, Aeronautics
Institute of Technology, 12228-900 São José dos Campos,
Brazil; ● orcid.org/0000-0002-6614-7625

Markus Oppel — Institute of Theoretical Chemistry, Faculty of
Chemistry, University of Vienna, 1090 Vienna, Austria;
● orcid.org/0000-0001-6363-2310

Carol A. Parish — Department of Chemistry, Gottwald Center
for the Sciences, University of Richmond, Richmond, Virginia
23173, United States; ● orcid.org/0000-0003-2878-3070

Jiri Pittner — J. Heyrovsky Institute of Physical Chemistry,
Academy of Sciences of the Czech Republic, 18223 Praha 8,
Czech Republic; ● orcid.org/0000-0003-0397-8410

Luan G. F. dos Santos — Department of Chemistry and
Biochemistry, Texas Tech University, Lubbock, Texas 79409,
United States; ● orcid.org/0000-0002-5979-148X

Armin Scrinzi — Department of Physics, Ludwig Maximilian
University, 80333 Munich, Germany

Mahesh K. Sit — J. Heyrovsky Institute of Physical Chemistry,
Academy of Sciences of the Czech Republic, 18223 Praha 8,
Czech Republic

Rene F. K. Spada — Departamento de Física, Instituto
Tecnológico de Aeronáutica, 12228-900 São José dos
Campos, São Paulo, Brazil

Mushir Thodika — National Renewable Energy Laboratory,
Golden, Colorado 80401, United States

Daniel C. A. Valente — Departamento de Física, Instituto
Tecnológico de Aeronáutica, 12228-900 São José dos
Campos, São Paulo, Brazil

Álvaro Vázquez-Mayagoitia — Computational Science
Division, Argonne National Laboratory, Lemont, Illinois
60439, United States; ● orcid.org/0000-0002-1415-6300

Elizete Ventura — Departamento de Química, CCEN,
Universidade Federal da Paraíba, 58059-900 João Pessoa,
Brazil

Julia Westermayr — Wilhelm-Ostwald-Institute for Physical
and Theoretical Chemistry, Leipzig University, 04103
Leipzig, Germany; ScaDS.AI (Center for Scalable Data
Analytics and Artificial Intelligence) Dresden/Leipzig, 04105
Leipzig, Germany; ● orcid.org/0000-0002-6531-0742

Aleksandr Zaichenko — Department of Chemistry, University
at Buffalo, State University of New York, Buffalo, New York
14260-3000, United States

Zhiyong Zhang — Stanford Research Computing, Stanford
University, Stanford, California 94305, United States

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jPCA.5c02047

Author Contributions

The manuscript was written through the contributions of all
authors. All authors have either contributed directly to the
development of the COLUMBUS package or contributed to its
progress through collaborations with the developer team,
submitting bug reports, feature requests, or developing new
innovative applications. All authors have given approval to the
final version of the manuscript.

Notes

The authors declare no competing financial interest.
## ACKNOWLEDGMENTS

F.P. acknowledges the use of the Lovelace HPC service at Loughborough University. H.L. acknowledges support from the U.S. National Science Foundation, grant No. CHE-2107923, and for computer time at the HPCC facilities of Texas Tech University. F.B.C.M., V.P.O., and J.C.V.C. acknowledge the financial support of the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (grants 307168/2022-0 and 164766/2021-9), the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) (grants 2018/13673-7 and 2022/16385-8), and the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) (grants 88887.877979/2023-00 and 88881.799245/2022-01). S.Mu. and M.B. acknowledge the funding provided by the European Research Council (ERC) under the Advanced Grant SubNano (Grant agreement No. 832237). They also acknowledge the Centre de Calcul Intensif d'Aix-Marseille for granting access to its high-performance computing resources. P.G.S. acknowledges support by the National Research, Development and Innovation Fund (NKFIF) of Hungary, Grant No. 142634. E.V., S.A.d.M., and R.L.R.A. thank the funding provided by the National Institute of Science and Technology on Molecular Sciences (INCT-CiMol), Grant CNPq 406804/2022-2. They also thank the Laboratório Nacional de Computação Científica—LNCC for the computational facilities of the Santos Dumont (SDumont) supercomputer cluster and the Brazilian agency FINEP for financial support. E.V. and S.A.d.M. also thank the financial support from CNPq (Grant Numbers 308371/2021-6 and 313398/2023-2). R.F.K.S. thanks the financial support from CNPq under grants 307846/2021-0 and 405562/2022-5. A.Z. and J.A. acknowledge financial support from the U.S. Department of Energy, Office of Basic Energy Sciences, Heavy Element Chemistry program, under grant DE-SC0001136, and thank the Center for Computational Research (CCR) at the University at Buffalo for providing computational resources. This research utilized resources of the Argonne Leadership Computing Facility, a U.S. Department of Energy (DOE) Office of Science user facility at Argonne National Laboratory A.V.M. was supported by the Office of Science, U.S. Department of Energy, under contract DE-AC02-06CH11357. R.S. was supported by the Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, U.S. Department of Energy, under Contract No. DE-AC02-06CH11357. J.P. and M.K.S. acknowledge the Czech Science Foundation Grant 23-06364S and the Ministry of Education, Youth, and Sports of the Czech Republic, project No. CZ.02.01.01/00/22_008/0004558, Co-funded by the European Union. Computer facilities provided by the program "Projects of Large Infrastructure for Research, Development, and Innovations" (no. LM2010005) and computer time provided by the IT4I supercomputing center (Project ID:90254) are also acknowledged. S.Ma. and T.M. acknowledge support from NSF CHE-2303111. This work in part was authored by the National Renewable Energy Laboratory for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Funding was provided by the Computational Chemical Sciences program within the Office of Basic Energy Sciences, U.S. Department of Energy. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article

for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

## REFERENCES

(1) Plasser, F.; Barbatti, M.; Aquino, A. J. A.; Lischka, H. Electronically Excited States and Photodynamics: A Continuing Challenge. *Theor. Chem. Acc.* **2012**, *131* (1), No. 1073.

(2) González, L.; Escudero, D.; Serrano-Andrés, L. Progress and Challenges in the Calculation of Electronic Excited States. *ChemPhysChem* **2012**, *13* (1), 28–51.

(3) Boggio-Pasqua, M.; Bearpark, M. J.; Klene, M.; Robb, M. A. A Computational Strategy for Geometry Optimization of Ionic and Covalent Excited States, Applied to Butadiene and Hexatriene. *J. Chem. Phys.* **2004**, *120* (17), 7849–7860.

(4) Borden, W. T.; Davidson, E. R. The Importance of Including Dynamic Electron Correlation in *Ab Initio* Calculations. *Acc. Chem. Res.* **1996**, *29* (2), 67–75.

(5) Roos, B. O.; Andersson, K.; Fülscher, M. P. Towards an Accurate Molecular Orbital Theory for Excited States: The Benzene Molecule. *Chem. Phys. Lett.* **1992**, *192* (1), 5–13.

(6) Grimme, S.; Hansen, A. A Practicable Real-Space Measure and Visualization of Static Electron-Correlation Effects. *Angew. Chem., Int. Ed.* **2015**, *54* (42), 12308–12313.

(7) Shavitt, I. The Graphical Unitary Group Approach and Its Application to Direct Configuration Interaction Calculations. In *The Unitary Group for the Evaluation of Electronic Energy Matrix Elements*; Hinze, J., Ed.; Lecture Notes in Chemistry; Springer: Berlin, 1981; pp 51–99.

(8) Shepard, R.; Shavitt, I.; Pitzer, R. M.; Comeau, D. C.; Pepper, M.; Lischka, H.; Szalay, P. G.; Ahlrichs, R.; Brown, F. B.; Zhao, J.-G. A Progress Report on the Status of the COLUMBUSMRCI Program System. *Int. J. Quantum Chem.* **1988**, *34* (S22), 149–165.

(9) Plasser, F.; Lischka, H. Multi-Reference Configuration Interaction. In *Quantum Chemistry and Dynamics of Excited States*; Wiley, 2020; Chapter 9, pp 277–297.

(10) Müller, T. Large-Scale Parallel Uncontracted Multireference-Averaged Quadratic Coupled Cluster: The Ground State of the Chromium Dimer Revisited. *J. Phys. Chem. A* **2009**, *113* (45), 12729–12740.

(11) Plasser, F.; Pašalić, H.; Gerzabek, M. H.; Libisch, F.; Reiter, R.; Burgdörfer, J.; Müller, T.; Shepard, R.; Lischka, H. The Multiradical Character of One- and Two-Dimensional Graphene Nanoribbons. *Angew. Chem., Int. Ed.* **2013**, *52* (9), 2581–2584.

(12) Lischka, H.; Shepard, R.; Müller, T.; Szalay, P. G.; Pitzer, R. M.; Aquino, A. J. A.; do Nascimento, M. M. A.; Barbatti, M.; Belcher, L. T.; Blaudeau, J.-P.; et al. The Generality of the GUGA MRCI Approach in COLUMBUS for Treating Complex Quantum Chemistry. *J. Chem. Phys.* **2020**, *152* (13), No. 134110.

(13) Yabushita, S.; Zhang, Z.; Pitzer, R. M. Spin-Orbit Configuration Interaction Using the Graphical Unitary Group Approach and Relativistic Core Potential and Spin-Orbit Operators. *J. Phys. Chem. A* **1999**, *103* (29), 5791–5800.

(14) Mai, S.; Müller, T.; Plasser, F.; Marquetand, P.; Lischka, H.; González, L. Perturbational Treatment of Spin-Orbit Coupling for Generally Applicable High-Level Multi-Reference Methods. *J. Chem. Phys.* **2014**, *141* (7), No. 074105.

(15) Langhoff, S. R.; Davidson, E. R. Configuration Interaction Calculations on the Nitrogen Molecule. *Int. J. Quantum Chem.* **1974**, *8* (1), 61–72.

(16) Szalay, P. G. Configuration Interaction: Corrections for Size-Consistency. In *Encyclopedia of Computational Chemistry*; Wiley, 2005.

(17) do Monte, S. A.; Spada, R. F. K.; Alves, R. L. R.; Belcher, L.; Shepard, R.; Lischka, H.; Plasser, F. Quantification of the Ionic Character of Multiconfigurational Wave Functions: The Qat Diagnostic. *J. Phys. Chem. A* **2023**, *127* (46), 9842–9852.
(18) Pople, J. A.; Seeger, R.; Krishnan, R. Variational Configuration Interaction Methods and Comparison with Perturbation Theory. *Int. J. Quantum Chem.* **2009**, *12*, 149–163.

(19) Davidson, E. R.; Silver, D. W. Size Consistency in the Dilute Helium Gas Electronic Structure. *Chem. Phys. Lett.* **1977**, *52* (3), 403–406.

(20) Szalay, P. G.; Bartlett, R. J. Multi-Reference Averaged Quadratic Coupled-Cluster Method: A Size-Extensive Modification of Multi-Reference CI. *Chem. Phys. Lett.* **1993**, *214* (5), 481–488.

(21) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; et al. The Dalton Quantum Chemistry Program System. *WIREs Comput. Mol. Sci.* **2014**, *4* (3), 269–284.

(22) Li Manni, G.; Fdez; Galván, I.; Alavi, A.; Aleotti, F.; Aquilante, F.; Autschbach, J.; Avagliano, D.; Baiardi, A.; Bao, J. J.; Battaglia, S.; et al. The OpenMolcas Web: A Community-Driven Approach to Advancing Computational Chemistry. *J. Chem. Theory Comput* **2023**, *19* (20), 6933–6991.

(23) Aprà, E.; Bylaska, E. J.; de Jong, W. A.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Valiev, M.; van Dam, H. J. J.; Alexeev, Y.; Anchell, J.; et al. NWChem: Past, Present, and Future. *J. Chem. Phys.* **2020**, *152* (18), No. 184102, DOI: 10.1063/5.0004997.

(24) Aquilante, F.; Autschbach, J.; Carlson, R. K.; Chibotaru, L. F.; Delcey, M. G.; De Vico, L.; Fdez; Galván, I.; Ferré, N.; Frutos, L. M.; Gagliardi, L.; et al. Molcas 8: New Capabilities for Multiconfigurational Quantum Chemical Calculations across the Periodic Table. *J. Comput. Chem.* **2016**, *37* (5), 506–541.

(25) COLUMBUS. an Ab Initio Electronic Structure Program. https://gitlab.com/columbus-program-system/columbus. (accessed June 09, 2025).

(26) Spada, R. F. K.; Franco, M. P.; Nieman, R.; Aquino, A. J. A.; Shepard, R.; Plasser, F.; Lischka, H. Spin-Density Calculation via the Graphical Unitary Group Approach. *Mol. Phys.* **2023**, *121* (11–12), No. e2091049.

(27) Shavitt, I. Graph Theoretical Concepts for the Unitary Group Approach to the Many-Electron Correlation Problem. *Int. J. Quantum Chem.* **2009**, *12* (S11), 131–148.

(28) Gidofalvi, G.; Shepard, R. The Evaluation of Spin-density Matrices within the Graphically Contracted Function Method. *Int. J. Quantum Chem.* **2009**, *109* (15), 3552–3563.

(29) Nakano, M.; Champagne, B. Theoretical Design of Open-Shell Singlet Molecular Systems for Nonlinear Optics. *J. Phys. Chem. Lett.* **2015**, *6* (16), 3236–3256.

(30) Messelberger, J.; Grünwald, A.; Pinter, P.; Hansmann, M. M.; Munz, D. Carbene Derived Diradicaloids – Building Blocks for Singlet Fission? *Chem. Sci.* **2018**, *9* (28), 6107–6117.

(31) Yu, C. P.; Chowdhury, R.; Fu, Y.; Ghosh, P.; Zeng, W.; Mustafa, T. B. E.; Grüne, J.; Walker, L. E.; Congrave, D. G.; Chua, X. W.; Murto, P.; Rao, A.; Sirringhaus, H.; Plasser, F.; Grey, C. P.; Friend, R. H.; Bronstein, H. Near-Infrared Luminescent Open-Shell π-Conjugated Systems with a Bright Lowest-Energy Zwitterionic Singlet Excited State. *Sci. Adv.* **2024**, *10* (30), No. eado3476, DOI: 10.1126/sciadv.ado3476.

(32) Rudebusch, G. E.; Zafra, J. L.; Jorner, K.; Fukuda, K.; Marshall, J. L.; Arrechea-Marcos, I.; Espejo, G. L.; Ponce Ortiz, R.; Gómez-García, C. J.; Zakharov, L. N.; Nakano, M.; Ottosson, H.; Casado, J.; Haley, M. M. Diindeno-Fusion of an Anthracene as a Design Strategy for Stable Organic Biradicals. *Nat. Chem.* **2016**, *8* (8), 753–759.

(33) Nguyen, T. D.; Ehrenfreund, E.; Vardeny, Z. V. Spin-Polarized Light-Emitting Diode Based on an Organic Bipolar Spin Valve. *Science* **2012**, *337* (6091), 204–209.

(34) Szwarc, M. Some Remarks on the CH2CH2Molecule. *Discuss. Faraday Soc.* **1947**, *2*, 46–49.

(35) Coulson, C. A.; Craig, D. P.; Maccoll, A.; Pullman, A. P.; Quinodimethane and Its Diradical. *Discuss. Faraday Soc.* **1947**, *2*, 36–38.

(36) Chagas, J. C. V.; Milanez, B. D.; Oliveira, V. P.; Pinheiro, M., Jr.; Ferrão, L. F. A.; Aquino, A. J. A.; Lischka, H.; Machado, F. B. C. A Multi-descriptor Analysis of Substituent Effects on the Structure and

Aromaticity of Benzene Derivatives: II-Conjugation versus Charge Effects. *J. Comput. Chem.* **2024**, *45* (12), 863–877.

(37) Matasović, L.; Bronstein, H.; Friend, R. H.; Plasser, F. Classification and Quantitative Characterisation of the Excited States of π-Conjugated Diradicals. *Faraday Discuss.* **2024**, *254*, 107–129.

(38) Schulten, K.; Ohmine, I.; Karplus, M. Correlation Effects in the Spectra of Polyenes. *J. Chem. Phys.* **1976**, *64* (11), 4422–4441.

(39) Kimber, P.; Plasser, F. Classification and Analysis of Molecular Excited States. In *Comprehensive Computational Chemistry*, Elsevier, 2024; pp 55–83.

(40) Plasser, F. TheoDORE: A Toolbox for a Detailed and Automated Analysis of Electronic Excited State Computations. *J. Chem. Phys.* **2020**, *152* (8), No. 084108.

(41) Véril, M.; Scemama, A.; Caffarel, M.; Lipparini, F.; Boggio-Pasqua, M.; Jacquemin, D.; Loos, P. QUESTDB: A Database of Highly Accurate Excitation Energies for the Electronic Structure Community. *WIREs Comput. Mol. Sci.* **2021**, *11* (5), No. e1517.

(42) Chen, X. K.; Kim, D.; Brédas, J. L. Thermally Activated Delayed Fluorescence (TADF) Path toward Efficient Electroluminescence in Purely Organic Materials: Molecular Level Insight. *Acc. Chem. Res.* **2018**, *51* (9), 2215–2224.

(43) Kimber, P.; Plasser, F. Energy Component Analysis for Electronically Excited States of Molecules: Why the Lowest Excited State Is Not Always the HOMO/LUMO Transition. *J. Chem. Theory Comput.* **2023**, *19* (8), 2340–2352.

(44) Kimber, P.; Plasser, F. Toward an Understanding of Electronic Excitation Energies beyond the Molecular Orbital Picture. *Phys. Chem. Chem. Phys.* **2020**, *22* (11), 6058–6080.

(45) Snyder, J. W.; Parrish, R. M.; Martínez, T. J. α-CASSCF: An Efficient, Empirical Correction for SA-CASSCF to Closely Approximate MS-CASPT2 Potential Energy Surfaces. *J. Phys. Chem. Lett.* **2017**, *8* (11), 2432–2437.

(46) Frutos, L. M.; Andruniów, T.; Santoro, F.; Ferré, N.; Olivucci, M. Tracking the Excited-State Time Evolution of the Visual Pigment with Multiconfigurational Quantum Chemistry. *Proc. Natl. Acad. Sci. U.S.A.* **2007**, *104* (19), 7764–7769.

(47) Carlson, R. K.; Truhlar, D. G.; Gagliardi, L. Multiconfiguration Pair-Density Functional Theory: A Fully Translated Gradient Approximation and Its Performance for Transition Metal Dimers and the Spectroscopy of Re2Cl82. *J. Chem. Theory Comput* **2015**, *11* (9), 4077–4085.

(48) Gdanitz, R. J.; Ahlrichs, R. The Averaged Coupled-Pair Functional (ACPF): A Size-Extensive Modification of MR CI(SD). *Chem. Phys. Lett.* **1988**, *143* (5), 413–420.

(49) Szalay, P. G.; Bartlett, R. J. Approximately Extensive Modifications of the Multireference Configuration Interaction Method: A Theoretical and Practical Analysis. *J. Chem. Phys.* **1995**, *103* (9), 3600–3612.

(50) Szalay, P. G.; Müller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration Self-Consistent Field and Multireference Configuration Interaction Methods and Applications. *Chem. Rev.* **2012**, *112* (1), 108–181.

(51) Lischka, H.; Nachtigallová, D.; Aquino, A. J. A.; Szalay, P. G.; Plasser, F.; Machado, F. B. C.; Barbatti, M. Multireference Approaches for Excited States of Molecules. *Chem. Rev.* **2018**, *118* (15), 7293–7361.

(52) Szalay, P. G. Towards State-Specific Formulation of Multi-reference Coupled-Cluster Theory: Coupled Electron Pair Approximations (CEPA) Leading to Multireference Configuration Interaction (MR-CI) Type Equations. In *Recent Advances in Coupled-Cluster Methods*, Bartlett, R. J., Ed.; World Scientific, 1997; pp 81–123.

(53) Szalay, P. G. Configuration Interaction: Corrections for Size-Consistency. In *Encyclopedia of Computational Chemistry*; Wiley, 1998.

(54) Szalay, P. G.; Müller, T.; Lischka, H. Excitation Energies and Transition Moments by the Multireference Averaged Quadratic Coupled Cluster (MR-AQCC) Method. *Phys. Chem. Chem. Phys.* **2000**, *2* (10), 2067–2073.
(55) Lischka, H.; Shepard, R.; Pitzer, R. M.; Shavitt, I.; Dallos, M.; Müller, T.; Szalay, P. G.; Seth, M.; Kedziora, G. S.; Yabushita, S.; Zhang, Z. High-Level Multireference Methods in the Quantum-Chemistry Program System COLUMBUS: Analytic MR-CISD and MR-AQCC Gradients and MR-AQCC-LRT for Excited States, GUGA Spin-Orbit CI and Parallel CI Density. *Phys. Chem. Chem. Phys.* **2001**, *3* (5), 664–673.

(56) Horn, S.; Plasser, F.; Müller, T.; Libisch, F.; Burgdörfer, J.; Lischka, H. A Comparison of Singlet and Triplet States for One- and Two-Dimensional Graphene Nanoribbons Using Multireference Theory. *Theor. Chem. Acc.* **2014**, *133* (8), No. 1511.

(57) Milanez, B. D.; Chagas, J. C. V.; Pinheiro, M., Jr; Aquino, A. J. A.; Lischka, H.; Machado, F. B. C. Effects on the Aromaticity and on the Biradicaloid Nature of Acenes by the Inclusion of a Cyclobutadiene Linkage. *Theor. Chem. Acc.* **2020**, *139* (7), No. 113.

(58) Müller, T.; Dallos, M.; Lischka, H.; Dubrovay, Z.; Szalay, P. G. A Systematic Theoretical Investigation of the Valence Excited States of the Diatomic Molecules B 2, C 2, N 2 and O 2. *Theor. Chem. Acc.: Theory, Comput., Model. (Theor. Chim. Acta)* **2001**, *105* (3), 227–243.

(59) Christiansen, O.; Jorgensen, P.; Hättig, C. Response Functions from Fourier Component Variational Perturbation Theory Applied to a Time-Averaged Quasienergy. *Int. J. Quantum Chem.* **1998**, *68* (1), 1–52.

(60) Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. *J. Chem. Phys.* **1989**, *90* (2), 1007–1023.

(61) Seritan, S.; Bannwarth, C.; Fales, B. S.; Hohenstein, E. G.; Isborn, C. M.; Kokkila-Schumacher, S. I. L.; Li, X.; Liu, F.; Luehr, N.; Snyder, J. W.; Song, C.; Titov, A. V.; Ufimtsev, I. S.; Wang, L.; Martínez, T. J. TeraChem: A Graphical Processing Unit-accelerated Electronic Structure Package for large-scale Ab Initio Molecular Dynamics. *WIREs Comput. Mol. Sci.* **2021**, *11* (2), No. e1494.

(62) NVIDIA Corporation. NVIDIA CUDA Fortran Programming Guide 2025. https://docs.nvidia.com/hpc-sdk/compilers/cuda-fortran-prog-guide/index.html; https://docs.nvidia.com/hpc-sdk/compilers/cuda-fortran-prog-guide/index.html. (accessed January 31, 2025).

(63) Wilkinson, K. A.; Sherwood, P.; Guest, M. F.; Naidoo, K. J. Acceleration of the GAMESS-UK Electronic Structure Package on Graphical Processing Units. *J. Comput. Chem.* **2011**, *32* (10), 2313–2318.

(64) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B. et al. Gaussian 16; Gaussian, Inc.: Wallingford, 2016.

(65) Sun, Q.; Zhang, X.; Banerjee, S.; Bao, P.; Barbry, M.; Blunt, N. S.; Bogdanov, N. A.; Booth, G. H.; Chen, J.; Cui, Z.-H.; Eriksen, J. J.; Gao, Y.; Guo, S.; Hermann, J.; Hermes, M. R.; Koh, K.; Koval, P.; Lehtola, S.; Li, Z.; Liu, J.; Mardirossian, N.; McClain, J. D.; Motta, M.; Mussard, B.; Pham, H. Q.; Pulkin, A.; Purwanto, W.; Robinson, P. J.; Ronca, E.; Sayfutyarova, E. R.; Scheurer, M.; Schurkus, H. F.; Smith, J. E. T.; Sun, C.; Sun, S.-N.; Upadhyay, S.; Wagner, L. K.; Wang, X.; White, A.; Whitfield, J. D.; Williamson, M. J.; Wouters, S.; Yang, J.; Yu, J. M.; Zhu, T.; Berkelbach, T. C.; Sharma, S.; Sokolov, A. Yu.; Chan, G. K.-L. Recent Developments in the PySCF Program Package. *J. Chem. Phys.* **2020**, *153* (2), No. 024109.

(66) Roos, B. O.; Siegbahn, P. E. M. The Direct Configuration Interaction Method from Molecular Integrals. In *Methods of Electronic Structure Theory*; Schaefer, H. F., III, Ed.; Springer: New York, 1977; pp 277–318.

(67) Paldus, J. Group Theoretical Approach to the Configuration Interaction and Perturbation Theory Calculations for Atomic and Molecular Systems. *J. Chem. Phys.* **1974**, *61* (12), 5321–5330.

(68) Ahlrichs, R.; Böhm, H.; Ehrhardt, C.; Scharf, P.; Schiffer, H.; Lischka, H.; Schindler, M. Implementation of an Electronic Structure Program System on the CYBER 205. *J. Comput. Chem.* **1985**, *6* (3), 200–208.

(69) Werner, H.-J.; Reinsch, E. The Self-consistent Electron Pairs Method for Multiconfiguration Reference State Functions. *J. Chem. Phys.* **1982**, *76* (6), 3144–3156.

(70) Davidson, E. R. The Iterative Calculation of a Few of the Lowest Eigenvalues and Corresponding Eigenvectors of Large Real-Symmetric Matrices. *J. Comput. Phys.* **1975**, *17* (1), 87–94.

(71) Shepard, R.; Shavitt, I.; Lischka, H. Reducing I/O Costs for the Eigenvalue Procedure in Large-scale Configuration Interaction Calculations. *J. Comput. Chem.* **2002**, *23* (11), 1121–1125.

(72) Arumainayagam, C. R.; Garrod, R. T.; Boyer, M. C.; Hay, A. K.; Bao, S. T.; Campbell, J. S.; Wang, J.; Nowak, C. M.; Arumainayagam, M. R.; Hodge, P. J. Extraterrestrial Prebiotic Molecules: Photo-chemistry vs. Radiation Chemistry of Interstellar Ices. *Chem. Soc. Rev.* **2019**, *48* (8), 2293–2314.

(73) Boyer, M. C.; Rivas, N.; Tran, A. A.; Verish, C. A.; Arumainayagam, C. R. The Role of Low-Energy (≤ 20 EV) Electrons in Astrochemistry. *Surf. Sci.* **2016**, *652*, 26–32.

(74) Swiderek, P. Fundamental Processes in Radiation Damage of DNA. *Angew. Chem., Int. Ed.* **2006**, *45* (25), 4056–4059.

(75) Alizadeh, E.; Orlando, T. M.; Sanche, L. Biomolecular Damage Induced by Ionizing Radiation: The Direct and Indirect Effects of Low-Energy Electrons on DNA. *Annu. Rev. Phys. Chem.* **2015**, *66* (1), 379–398.

(76) Arumainayagam, C. R.; Lee, H.-L.; Nelson, R. B.; Haines, D. R.; Gunawardane, R. P. Low-Energy Electron-Induced Reactions in Condensed Matter. *Surf. Sci. Rep* **2010**, *65* (1), 1–44.

(77) Herbert, J. M. The Quantum Chemistry of Loosely-Bound Electrons. In *Reviews in Computational Chemistry*; Parrill, A. L.; Lipkowitz, K. B., Eds.; Wiley, 2015;Chapter 8, Vol. 28, pp 391–517.

(78) Kawarai, Y.; Weber, Th.; Azuma, Y.; Winstead, C.; McKoy, V.; Belkacem, A.; Slaughter, D. S. Dynamics of the Dissociating Uracil Anion Following Resonant Electron Attachment. *J. Phys. Chem. Lett.* **2014**, *5* (21), 3854–3858.

(79) Klaiman, S.; Cederbaum, L. S. Barrierless Single-Electron-Induced cis-rans Isomerization. *Angew. Chem.* **2015**, *127* (36), 10616–10619.

(80) Fennimore, M. A.; Karsili, T. N. V.; Matsika, S. Mechanisms of H and CO Loss from the Uracil Nucleobase Following Low Energy Electron Irradiation. *Phys. Chem. Chem. Phys.* **2017**, *19* (26), 17233–17241.

(81) Fennimore, M. A.; Matsika, S. Electronic Resonances of Nucleobases Using Stabilization Methods. *J. Phys. Chem. A* **2018**, *122* (16), 4048–4057.

(82) Loupas, A.; Gorfinkiel, J. D. Shape and Core-Excited Resonances in Electron Scattering from Alanine. *J. Chem. Phys.* **2019**, *150* (6), No. 064307.

(83) Thodika, M.; Mackouse, N.; Matsika, S. Description of Two-Particle One-Hole Electronic Resonances Using Orbital Stabilization Methods. *J. Phys. Chem. A* **2020**, *124* (43), 9011–9020.

(84) Gayvert, J. R. An Open-Source Program for Studying Resonances in Molecules 2021 https://github.com/gayverjr/opencap. (accessed June 08, 2025).

(85) Moiseyev, N. *Non-Hermitian Quantum Mechanics*; Cambridge University Press, 2011.

(86) Riss, U. V.; Meyer, H.-D. Calculation of Resonance Energies and Widths Using the Complex Absorbing Potential Method. *J. Phys. B: At., Mol. Opt. Phys.* **1993**, *26* (23), 4503–4535.

(87) Jagau, T.-C.; Zuev, D.; Bravaya, K. B.; Epifanovsky, E.; Krylov, A. I. A Fresh Look at Resonances and Complex Absorbing Potentials: Density Matrix-Based Approach. *J. Phys. Chem. Lett.* **2014**, *5* (2), 310–315.

(88) Sommerfeld, T.; Santra, R. Efficient Method to Perform CAP/CI Calculations for Temporary Anions. *Int. J. Quantum Chem.* **2001**, *82* (5), 218–226.

(89) Santra, R.; Cederbaum, L. S. An Efficient Combination of Computational Techniques for Investigating Electronic Resonance States in Molecules. *J. Chem. Phys.* **2001**, *115* (15), 6853–6861.

(90) Thodika, M.; Matsika, S. Projected Complex Absorbing Potential Multireference Configuration Interaction Approach for Shape and Feshbach Resonances. *J. Chem. Theory Comput* **2022**, *18* (6), 3377–3390.
(91) Thodika, M. *Development and Benchmarking of Hermitian and Non-Hermitian Methods for Negative Ion Resonances*; Temple University, 2022.

(92) Majety, V. P.; Zielinski, A.; Scrinzi, A. Photoionization of Few Electron Systems: A Hybrid Coupled Channels Approach. *New J. Phys.* **2015**, *17* (6), No. 063002.

(93) Chundayil, H.; Majety, V. P.; Scrinzi, A. The Hybrid Anti-Symmetrized Coupled Channels Method (HaCC) for the TRecX Code. *Comput. Phys. Commun.* **2024**, *303*, No. 109279.

(94) Majety, V. P.; Scrinzi, A. Dynamic Exchange in the Strong Field Ionization of Molecules. *Phys. Rev. Lett.* **2015**, *115* (10), No. 103002.

(95) Tao, L.; Scrinzi, A. Photo-Electron Momentum Spectra from Minimal Volumes: The Time-Dependent Surface Flux Method. *New J. Phys.* **2012**, *14* (1), No. 013021.

(96) Majety, V. P.; Scrinzi, A. Multielectron Effects in Strong-Field Ionization of CO2: Impact on Differential Photoelectron Spectra. *Phys. Rev. A* **2017**, *96* (5), No. 053421.

(97) Geim, A. K.; Novoselov, K. S. The Rise of Graphene. *Nat. Mater.* **2007**, *6* (3), 183–191.

(98) Nguyen, B. H.; Nguyen, V. H. Promising Applications of Graphene and Graphene-Based Nanostructures. *Adv. Nat. Sci.: Nanosci. Nanotechnol.* **2016**, *7* (2), No. 023002.

(99) Gu, Y.; Qiu, Z.; Müllen, K. Nanographenes and Graphene Nanoribbons as Multitalents of Present and Future Materials Science. *J. Am. Chem. Soc.* **2022**, *144* (26), 11499–11524.

(100) Nieman, R.; Silva, N. J.; Aquino, A. J. A.; Haley, M. M.; Lischka, H. Interplay of Biradicaloid Character and Singlet/Triplet Energy. Splitting for *Cis* / *Trans* -Diindenoacenes and Related Benzothiophene-Capped Oligomers as Revealed by Extended Multi-reference Calculations. *J. Org. Chem.* **2020**, *85* (5), 3664–3675.

(101) Kubo, T. Phenalenyl-Based Open-Shell Polycyclic Aromatic Hydrocarbons. *Chem. Rec.* **2015**, *15* (1), 218–232.

(102) Sánchez-Grande, A.; Urgel, J. I.; Cahlík, A.; Santos, J.; Edalatmanesh, S.; Rodríguez-Sánchez, E.; Lauwaet, K.; Mutombo, P.; Nachtigallová, D.; Nieman, R.; Lischka, H.; de la Torre, B.; Miranda, R.; Gröning, O.; Martín, N.; Jelínek, P.; Ecija, D. Diradical Organic One-Dimensional Polymers Synthesized on a Metallic Surface. *Angew. Chem., Int. Ed.* **2020**, *59* (40), 17594–17599.

(103) Bettinger, H. F. Electronic Structure of Higher Acenes and Polyacene: The Perspective Developed by Theoretical Analyses. *Pure Appl. Chem.* **2010**, *82* (4), 905–915.

(104) Ahmed, J.; Mandal, S. K. Phenalenyl Radical: Smallest Polycyclic Odd Alternant Hydrocarbon Present in the Graphene Sheet. *Chem. Rev.* **2022**, *122* (13), 11369–11431.

(105) Kubo, T. Closed-Shell and Open-Shell Dual Nature of Singlet Diradical Compounds. *Pure Appl. Chem.* **2023**, *95* (4), 363–375.

(106) Takatsuka, K.; Fueno, T.; Yamaguchi, K. Distribution of Odd Electrons in Ground-State Molecules. *Theor. Chim. Acta* **1978**, *48* (3), 175–183.

(107) Head-Gordon, M. Characterizing Unpaired Electrons from the One-Particle Density Matrix. *Chem. Phys. Lett.* **2003**, *372* (3–4), 508–511.

(108) Plasser, F.; Wormit, M.; Dreuw, A. New Tools for the Systematic Analysis and Visualization of Electronic Excitations. I. Formalism. *J. Chem. Phys.* **2014**, *141* (2), No. 024106.

(109) Bauer, C. A.; Hansen, A.; Grimme, S. The Fractional Occupation Number Weighted Density as a Versatile Analysis Tool for Molecules with a Complicated Electronic Structure. *Chem. Eur. J.* **2017**, *23* (25), 6150–6164.

(110) Nieman, R.; Carvalho, J. R.; Jayee, B.; Hansen, A.; Aquino, A. J. A.; Kertesz, M.; Lischka, H. Polyradical Character Assessment Using Multireference Calculations and Comparison with Density-Functional Derived Fractional Occupation Number Weighted Density Analysis. *Phys. Chem. Chem. Phys.* **2023**, *25* (40), 27380–27393.

(111) Carvalho, J. R.; Nieman, R.; Kertesz, M.; Aquino, A. J. A.; Hansen, A.; Lischka, H. Multireference Calculations on Bond Dissociation and Biradical Polycyclic Aromatic Hydrocarbons as Guidance for Fractional Occupation Number Weighted Density

Analysis in DFT Calculations. *Theor. Chem. Acc.* **2024**, *143* (10), No. 69.

(112) Staroverov, V. N.; Davidson, E. R. Distribution of Effectively Unpaired Electrons. *Chem. Phys. Lett.* **2000**, *330* (1–2), 161–168.

(113) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids. *Phys. Rev. Lett.* **2003**, *91* (14), 3–6.

(114) 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* (18), 3297–3305.

(115) Tsuneda, T.; Hirao, K. Long-range Correction for Density Functional Theory. *WIREs Comput. Mol. Sci.* **2014**, *4* (4), 375–390.

(116) Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-XTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions. *J. Chem. Theory Comput* **2019**, *15* (3), 1652–1671.

(117) 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* (2), No. e1493.

(118) Bhatt, M. D.; Kim, H.; Kim, G. Various Defects in Graphene: A Review. *RSC Adv.* **2022**, *12* (33), 21520–21547.

(119) Machado, F. B. C.; Aquino, A. J. A.; Lischka, H. The Diverse Manifold of Electronic States Generated by a Single Carbon Defect in a Graphene Sheet: Multireference Calculations Using a Pyrene Defect Model. *ChemPhysChem* **2014**, *15* (15), 3334–3341.

(120) Pinheiro, M.; Cardoso, D. V. V.; Aquino, A. J. A.; Machado, F. B. C.; Lischka, H. The Characterization of Electronic Defect States of Single and Double Carbon Vacancies in Graphene Sheets Using Molecular Density Functional Theory. *Mol. Phys.* **2019**, *117* (9–12), 1519–1531.

(121) Nieman, R.; Das, A.; Aquino, A. J. A.; Amorim, R. G.; Machado, F. B. C.; Lischka, H. Single and Double Carbon Vacancies in Pyrene as First Models for Graphene Defects: A Survey of the Chemical Reactivity toward Hydrogen. *Chem. Phys.* **2017**, *482*, 346–354.

(122) Nieman, R.; Oliveira, V. P.; Jayee, B.; Adelia, A. J. A.; Machado, F. B. C.; Lischka, H. High-Level Multireference Investigations on the Electronic States in Single-Vacancy (SV) Graphene Defects Using a Pyrene-SV Model. *J. Phys. Chem. A* **2023**, *127* (40), 8287–8296.

(123) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. *Theor. Chim. Acta* **1973**, *28* (3), 213–222.

(124) Shu, C.; Yang, Z.; Rajca, A. From Stable Radicals to Thermally Robust High-Spin Diradicals and Triradicals. *Chem. Rev.* **2023**, *123* (20), 11954–12003.

(125) Pozo, I.; Bogani, L. A Perspective on Radicaloid Conjugated Polycyclic Hydrocarbons. *Trends Chem.* **2024**, *6* (10), 581–595.

(126) Zeng, W.; Wu, J. Open-Shell Graphene Fragments. *Chem* **2021**, *7* (2), 358–386.

(127) Yu, H.; Jing, Y.; Heine, T. Physics and Chemistry of Two-Dimensional Triangulene-Based Lattices. *Acc. Chem. Res.* **2025**, *58* (1), 61–72.

(128) Mishra, S.; Beyer, D.; Eimre, K.; Liu, J.; Berger, R.; Gröning, O.; Pignedoli, C. A.; Müllen, K.; Fasel, R.; Feng, X.; Ruffieux, P. Synthesis and Characterization of π-Extended Triangulene. *J. Am. Chem. Soc.* **2019**, *141* (27), 10621–10625.

(129) Sanz-Rodrigo, J.; Ricci, G.; Olivier, Y.; Sancho-García, J. C. Negative Singlet–Triplet Excitation Energy Gap in Triangle-Shaped Molecular Emitters for Efficient Triplet Harvesting. *J. Phys. Chem. A* **2021**, *125* (2), 513–522.

(130) Martínez-Carracedo, G.; Oroszlány, L.; García-Fuente, A.; Szunyogh, L.; Ferrer, J. Electrically Driven Singlet-Triplet Transition in Triangulene Spin-1 Chains. *Phys. Rev. B* **2023**, *107* (3), No. 035432.
(131) Su, J.; Telychko, M.; Song, S.; Lu, J. Triangulenes: From Precursor Design to On-Surface Synthesis and Characterization. *Angew. Chem., Int. Ed.* **2020**, *59* (20), 7658–7668.

(132) Pavlíček, N.; Mistry, A.; Majzik, Z.; Moll, N.; Meyer, G.; Fox, D. J.; Gross, L. Synthesis and Characterization of Triangulene. *Nat. Nanotechnol.* **2017**, *12* (4), 308–311.

(133) Su, J.; Fan, W.; Mutombo, P.; Peng, X.; Song, S.; Ondráček, M.; Golub, P.; Brabec, J.; Veis, L.; Telychko, M.; Jelínek, P.; Wu, J.; Lu, J. On-Surface Synthesis and Characterization of [7]Triangulene Quantum Ring. *Nano Lett.* **2021**, *21* (1), 861–867.

(134) Mishra, S.; Beyer, D.; Eimre, K.; Ortiz, R.; Fernández-Rossier, J.; Berger, R.; Gröning, O.; Pignedoli, C. A.; Fasel, R.; Feng, X.; Ruffieux, P. Collective All-Carbon Magnetism in Triangulene Dimers. *Angew. Chem., Int. Ed.* **2020**, *59* (29), 12041–12047.

(135) Yu, H.; Heine, T. Magnetic Coupling Control in Triangulene Dimers. *J. Am. Chem. Soc.* **2023**, *145* (35), 19303–19311.

(136) Yu, H.; Sun, J.; Heine, T. Predicting Magnetic Coupling and Spin-Polarization Energy in Triangulene Analogues. *J. Chem. Theory Comput.* **2023**, *19* (12), 3486–3497.

(137) Saleem, Y.; Steenbock, T.; Alhadi, E. R. J.; Pasek, W.; Bester, G.; Potasz, P. Superexchange Mechanism in Coupled Triangulenes Forming Spin-1 Chains. *Nano Lett.* **2024**, *24* (24), 7417–7423.

(138) Weng, T.; Xu, Z.; Li, K.; Guo, Y.; Chen, X.; Li, Z.; Sun, Z. 1,1'-Biolympicenyl: A Stable Non-Kekulé Diradical with a Small Singlet and Triplet Energy Gap. *J. Am. Chem. Soc.* **2024**, *146* (38), 26454–26465.

(139) Ortiz, R.; Catarina, G.; Fernández-Rossier, J. Theory of Triangulene Two-Dimensional Crystals. *2d Mater.* **2023**, *10* (1), No. 015015.

(140) Ovchinnikov, A. A. Multiplicity of the Ground State of Large Alternant Organic Molecules with Conjugated Bonds - (Do Organic Ferromagnetics Exist?). *Theor Chim Acta* **1978**, *47* (4), 297–304.

(141) Lieb, E. H. Two Theorems on the Hubbard Model. *Phys. Rev. Lett.* **1989**, *62* (10), 1201–1204.

(142) Itkis, M. E.; Chi, X.; Cordes, A. W.; Haddon, R. C. Magneto-Opto-Electronic Bistability in a Phenalenyl-Based Neutral Radical. *Science* **2002**, *296* (5572), 1443–1445.

(143) Yang, Y.; Blacque, O.; Sato, S.; Juríček, M. Cycloparaphenylene–Phenalenyl Radical and Its Dimeric Double Nano-hoop. *Angew. Chem., Int. Ed.* **2021**, *60* (24), 13529–13535.

(144) Pal, S. K.; Itkis, M. E.; Tham, F. S.; Reed, R. W.; Oakley, R. T.; Haddon, R. C. Resonating Valence-Bond Ground State in a Phenalenyl-Based Neutral Radical Conductor. *Science* **2005**, *309* (5732), 281–284.

(145) Pariyar, A.; Vijaykumar, G.; Bhunia, M.; Dey, S. Kr.; Singh, S. K.; Kurungot, S.; Mandal, S. K. Switching Closed-Shell to Open-Shell Phenalenyl: Toward Designing Electroactive Materials. *J. Am. Chem. Soc.* **2015**, *137* (18), 5955–5960.

(146) Morita, Y.; Aoki, T.; Fukui, K.; Nakazawa, S.; Tamaki, K.; Suzuki, S.; Fuyuhiro, A.; Yamamoto, K.; Sato, K.; Shiomi, D.; Naito, A.; Takui, T.; Nakasuji, K. A New Trend in Phenalenyl Chemistry: A Persistent Neutral Radical, 2,5,8-Tri-Tert-Butyl-1,3-Diazaphenalenyl, and the Excited Triplet State of the Gablesyn-Dimer in the Crystal of Column Motif. *Angew. Chem., Int. Ed.* **2002**, *41* (10), 1793–1796.

(147) Das, A.; Müller, T.; Plasser, F.; Lischka, H. Polyradical Character of Triangular Non-Kekulé Structures, Zethrenes, p-Quinodimethane-Linked Bisphenalenyl, and the Clar Goblet in Comparison: An Extended Multireference Study. *J. Phys. Chem. A* **2016**, *120* (9), 1625–1636.

(148) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. *J. Chem. Phys.* **1972**, *56* (5), 2257–2261.

(149) Grimme, S. Semiempirical GGA-type Density Functional Constructed with a Long-range Dispersion Correction. *J. Comput. Chem.* **2006**, *27* (15), 1787–1799.

(150) Sirianni, D. A.; Song, X.; Wairegi, S.; Wang, E. B.; Mendoza-Gomez, S. A.; Luxon, A.; Zimmerley, M.; Nussdorf, A.; Filatov, M.; Hoffmann, R.; Parish, C. A. Variations on the Bergman Cyclization

Theme: Electrocyclizations of Ionic Penta-, Hepta-, and Octadiynes. *J. Am. Chem. Soc.* **2023**, *145* (39), 21408–21418.

(151) Jones, R. R.; Bergman, R. G. P-Benzyne. Generation as an Intermediate in a Thermal Isomerization Reaction and Trapping Evidence for the 1,4-Benzenediyl Structure. *J. Am. Chem. Soc.* **1972**, *94* (2), 660–661.

(152) Church, D. F.; Pryor, W. A. Free-Radical Chemistry of Cigarette Smoke and Its Toxicological Implications. *Environ. Health Perspect.* **1985**, *64*, 111–126.

(153) Russell, K. E.; Tobolsky, A. V. Diradicals in Solution: Role in Polymerization. *J. Am. Chem. Soc.* **1954**, *76* (2), 395–399.

(154) Slipchenko, L. V.; Krylov, A. I. Singlet-Triplet Gaps in Diradicals by the Spin-Flip Approach: A Benchmark Study. *J. Chem. Phys.* **2002**, *117* (10), 4694–4708.

(155) Wang, E. B.; Parish, C. A.; Lischka, H. An Extended Multireference Study of the Electronic States of *Para*-Benzyne. *J. Chem. Phys.* **2008**, *129* (4), No. 044306.

(156) Vu, K.; Pandian, J.; Zhang, B.; Annas, C.; Parker, A. J.; Mancini, J. S.; Wang, E. B.; Saldana-Greco, D.; Nelson, E. S.; Springsted, G.; Lischka, H.; Plasser, F.; Parish, C. A. Multireference Averaged Quadratic Coupled Cluster (MR-AQCC) Study of the Geometries and Energies for *Ortho*-, *Meta* - and *Para*-Benzyne. *J. Phys. Chem. A* **2024**, *128* (37), 7816–7829.

(157) Scott, T.; Nieman, R.; Luxon, A.; Zhang, B.; Lischka, H.; Gagliardi, L.; Parish, C. A. A Multireference Ab Initio Study of the Diradical Isomers of Pyrazine. *J. Phys. Chem. A* **2019**, *123* (10), 2049–2057.

(158) Pandian, J.; Vu, K.; Muya, J. T.; Parker, A.; Ancajas, C. M. F.; Saldana-Greco, D.; Yewer, T.; Parish, C. A. Highly Correlated, Multireference Study of the Lowest Lying Singlet and Triplet States of the Four Thiophene Diradicals. *J. Comput. Chem.* **2025**, *46* (3), No. e70044.

(159) Okamoto, K.; Kitagawa, T.; Takeuchi, K.; Komatsu, K.; Takahashi, K. Isolation of a Hydrocarbon Salt and Preparation of a Hydrocarbon Which Heterolyses to a Carbocation and a Carbanion. *J. Chem. Soc. Chem. Commun.* **1985**, No. 3, 173–174.

(160) Okamoto, K.; Kitagawa, T.; Takeuchi, K.; Komatsu, K.; Miyabo, A. A Hydrocarbon Existing Uniquely in Solution: A Heterolytically Dissociative Hydrocarbon That Produces the Corresponding Hydrocarbon Salt by Crystallization. *J. Chem. Soc. Chem. Commun.* **1988**, No. 14, 923–924.

(161) Kitagawa, T.; Tanaka, T.; Murakita, H.; Nishikawa, A.; Takeuchi, K. Reaction of Cyclopropenylium Ions with the Tert-Butyl-C60 Anion: Carbocation–Carbanion Coordination vs Salt Formation. *Tetrahedron* **2001**, *57* (17), 3537–3547.

(162) Tanaka, T.; Kitagawa, T.; Komatsu, K.; Takeuchi, K. Synthesis of a Hydrocarbon Salt Having a Fullerene Framework. *J. Am. Chem. Soc.* **1997**, *119* (39), 9313–9314.

(163) Ventura, E.; Alves, R. L. R.; do Monte, S. A. The Kinetics of Three Coupled Irreversible Elementary Reactions: Two Parallel Mixed Second Order Reactions Followed by a First Order Reaction. *J. Math. Chem.* **2024**, *62* (4), 922–935.

(164) Ventura, E.; Rodrigues, G. P.; Leitão, E. F. V.; do Monte, S. A. Theoretical Study of an Authentic Hydrocarbon Ion Pair. *ACS Omega* **2024**, *9* (32), 34981–34989.

(165) Alves, R. L. R.; Leitão, E. F. V.; Ventura, E.; do Monte, S. A. A Genuine Hydrocarbon Ion Pair More Stable Than Its Covalent Counterpart. A Computational Study. *J. Comput. Chem.* **2025**, *46* (8), No. e70079, DOI: 10.1002/jcc.70079.

(166) Kolomnikova, G. D.; Parnes, Z. N. Advances in the Chemistry of the Tropylium Ion. *Russ. Chem. Rev.* **1967**, *36* (10), 735–753.

(167) Koenig, T.; Chang, J. C. Helium(I) Photoelectron Spectrum of Tropyl Radical. *J. Am. Chem. Soc.* **1978**, *100* (7), 2240–2242.

(168) McDonald, R. N.; Bianchina, E. J.; Tung, C. C. Electron Photodetachment of Cyclopentadienylidene Anion Radical in a Flowing Afterglow Apparatus: EA and.DELTA.Hf.Degree. of Cyclopentadienylidene. *J. Am. Chem. Soc.* **1991**, *113* (19), 7115–7121.