* (118) Kapusta, J.I.; Kharzeev, D.; McLerran, L.D. The Return of the prodigal Goldstone boson. _Phys. Rev. D_ 1996, _53_ , 5028.
  * (119) Csorgo, T.; Vertesi, R.; Sziklai, J. Indirect observation of an in-medium $\eta$’ mass reduction in $\sqrt{s_{NN}}=200$ GeV Au+Au collisions. _Phys. Rev. Lett._ 2010, _105_ , 182301.
  * (120) Kotov, A.Y.; Lombardo, M.P.; Trunin, A.M. Fate of the $\eta?$ in the quark gluon plasma. _Phys. Lett. B_ 2019, _794_ , 83.
  * (121) Shuryak, E.V. Which chiral symmetry is restored in hot QCD? _Comments Nucl. Part. Phys._ 1994, _21_ , 235.
  * (122) Cohen, T.D. The High temperature phase of QCD and U(1)-A symmetry. _Phys. Rev. D_ 1996, _54_ , R1867.
  * (123) Lee, S.H.; Hatsuda, T. U-a(1) symmetry restoration in QCD with N(f) flavors. _Phys. Rev. D_ 1996, _54_ , R1871.
  * (124) Meggiolaro, E.; Morda, A. Remarks on the $U(1)$ axial symmetry and the chiral transition in QCD at finite temperature. _Phys. Rev. D_ 2013, _88_ , 096010.
  * (125) Heller, M.; Mitter, M. Pion and $\eta$-meson mass splitting at the two-flavor chiral crossover. _Phys. Rev. D_ 2016, _94_ , 074002.
  * (126) Ishii, M.; Yonemura, K.; Takahashi, J.; Kouno, H.; Yahiro, M. Determination of $U(1)_{\rm A}$ restoration from pion and $a_{0}$-meson screening masses: Toward the chiral regime. _Phys. Rev. D_ 2016, _93_ , 016002.
  * (127) Glozman, L. Chiralspin Symmetry and Its Implications for QCD. _Universe_ 2019, _5_ , 38.
  * (128) Gómez Nicola, A.; de Elvira, J.R. Pseudoscalar susceptibilities and quark condensates: Chiral restoration and lattice screening masses. _JHEP_ 2016, _1603_ , 186.
  * (129) Azcoiti, V. Topology in the SU(Nf) chiral symmetry restored phase of unquenched QCD and axion cosmology. _Phys. Rev. D_ 2016, _94_ , 094505.
  * (130) Gómez Nicola, A.; de Elvira, J.R. Patterns and partners for chiral symmetry restoration. _Phys. Rev. D_ 2018, _97_ , 074016.
  * (131) Gómez Nicola, A.; Elvira, J.R.D. Chiral and $U(1)_{A}$ restoration for the scalar and pseudoscalar meson nonets. _Phys. Rev. D_ 2018, _98_ , 014020.
  * (132) Buchoff, M.I.; Cheng, M.; Christ, N.H.; Ding, H.-T.; Jung, C.; Karsch, F.; Lin, Z.; Mawhinney, R.D.; Mukherjee, S.; Petreczky, P.; et al. QCD chiral transition, U(1)A symmetry and the Dirac spectrum using domain wall fermions. _Phys. Rev. D_ 2014, _89_ , 054514.
  * (133) Aoki, S.; Fukaya, H.; Taniguchi, Y. Chiral symmetry restoration, eigenvalue density of Dirac operator and axial U(1) anomaly at finite temperature. _Phys. Rev. D_ 2012, _86_ , 114512.
  * (134) Cossu, G.; Aoki, S.; Fukaya, H.; Hashimoto, S.; Kaneko, T.; Matsufuru, H.; Noaki, J.I. Finite temperature study of the axial U(1) symmetry on the lattice with overlap fermion formulation. _Phys. Rev. D_ 2013, _87_ , 114514. Erratum in 2013, _88_ , 019901.
  * (135) Tomiya, A.; Cossu, G.; Aoki, S.; Fukaya, H.; Hashimoto, S.; Kaneko, T.; Noaki, J. Evidence of effective axial U(1) symmetry restoration at high temperature QCD. _Phys. Rev. D_ 2017, _96_ , 034509.
  * (136) Brandt, B.B.; Francis, A.; Meyer, H.B.; Philipsen, O.; Robaina, D.; Wittig, H. On the strength of the $U(1)_{A}$ anomaly at the chiral phase transition in $N_{f}=2$ QCD. _JHEP_ 2016, _1612_ , 158.
  * (137) Brandt, B.B.; Cè, M.; Francis, A.; Harris, T.; Meyer, H.B.; Wittig, H.; Philipsen, O. Testing the strength of the $\text{U}_{A}(1)$ anomaly at the chiral phase transition in two-flavour QCD. _arXiv_ 2019, arXiv:1904.02384.
  * (138) Lee, J.W.; Lucini, B.; Piai, M. Symmetry restoration at high-temperature in two-color and two-flavor lattice gauge theories. _JHEP_ 2017, _1704_ , 036.
  * (139) Chandrasekharan, S.; Li, A. Anomaly and a QCD-like phase diagram with massive bosonic baryons. _JHEP_ 2010, _1012_ , 021.
  * (140) Aarts, G.; Allton, C.; Hands, S.; Jäger, B.; Praki, C.; Skullerud, J.I. ‘Nucleons and parity doubling across the deconfinement transition. _Phys. Rev. D_ 2015, _92_ , 014503.
  * (141) Aarts, G.; Allton, C.; Boni, D.D.; Hands, S.; Jäger, B.; Praki, C.; Skullerud, J.I. Light baryons below and above the deconfinement transition: Medium effects and parity doubling. _JHEP_ 2017, _1706_ , 034.
  * (142) Bochicchio, M.; Maiani, L.; Martinelli, G.; Rossi, G.C.; Testa, M. Chiral Symmetry on the Lattice with Wilson Fermions. _Nucl. Phys. B_ 1985, _262_ , 331.
  * (143) Boucaud, P.; Leroy, J.P.; Yaouanc, A.L.; Micheli, J.; Pene, O.; Rodriguez-Quintero, J. Quark pseudoscalar vertex and quark mass function with clover fermions: Spontaneous symmetry breaking, OPE, symmetry restoration at small volume. _Phys. Rev. D_ 2010, _81_ , 094504.
  * (144) Karsch, F.; Laermann, E. Thermodynamics and in medium hadron properties from lattice QCD. _arXiv_ 2003, arXiv:hep-lat/0305025
  * (145) Ishii, M.; Kouno, H.; Yahiro, M. Model prediction for temperature dependence of meson pole masses from lattice QCD results on meson screening masses. _Phys. Rev. D_ 2017, _95_ , 114022.
  * (146) Cheng, M.; Datta, S.; Francis, A.; van der Heide, J.; Jung, C.; Kaczmarek, O.; Karsch, F.; Laermann, E.; Mawhinney, R.D.; Miao, C.; et al. Meson screening masses from lattice QCD with two light and the strange quark. _Eur. Phys. J. C_ 2011, _71_ , 1564.
  * (147) Cheng, M.; Datta, S.; Francis, A.; van der Heide, J.; Jung, C.; Kaczmarek, O.; Karsch, F.; Laermann, E.; Mawhinney, R.D.; Miao, C.; et al. The QCD equation of state with almost physical quark masses. _Phys. Rev. D_ 2008, _77_ , 014511.
  * (148) Maezawa, Y.; Bazavov, A.; Karsch, F.; Petreczky, P.; Mukherjee, S. Meson screening masses at finite temperature with Highly Improved Staggered Quarks. _arXiv_ 2013, arXiv:1312.4375
  * (149) Gómez Nicola, A.; Elvira, J.R.D.; Vioque-Rodríguez, A. The QCD topological charge and its thermal dependence: The role of the $\eta^{\prime}$. _JHEP_ 2019, _11_ , 086.
  * (150) Guo, X.K.; Guo, Z.H.; Oller, J.A.; Sanz-Cillero, J.J. Scrutinizing the $\eta$-$\eta^{\prime}$ mixing, masses and pseudoscalar decay constants in the framework of U(3) chiral effective field theory. _JHEP_ 2015, _1506_ , 175.
  * (151) di Cortona, G.G.; Hardy, E.; Vega, J.P.; Villadoro, G. The QCD axion, precisely. _JHEP_ 2016, _1601_ , 034.
  * (152) Bernard, V.; Descotes-Genon, S.; Toucas, G. Topological susceptibility on the lattice and the three-flavour quark condensate. _JHEP_ 2012, _1206_ , 051.
  * (153) Aoki, S.; Aoki, Y.; Bernard, C.; Blum, T.; Colangelo, G.; Della Morte, M.; Durr, S.; El-Khadra, A.X.; Fukaya, H.; Horsley, R.; et al. Review of lattice results concerning low-energy particle physics. _Eur. Phys. J. C_ 2017, _77_ , 112.
  * (154) Burger, F.; Ilgenfritz, E.M.; Lombardo, M.P.; Trunin, A. Chiral observables and topology in hot QCD with two families of quarks. _Phys. Rev. D_ 2018, _98_ , 094501.
  * (155) Leutwyler, H.; Smilga, A.V. Spectrum of Dirac operator and role of winding number in QCD. _Phys. Rev. D_ 1992, _46_ , 5607.
  * (156) Veneziano, G. U(1) Without Instantons. _Nucl. Phys. B_ 1979, _159_ , 213.
  * (157) Mao, Y-Y.; Chiu, T-W. Topological Susceptibility to the One-Loop Order in Chiral Perturbation Theory. _Phys. Rev. D_ 2009, _80_ , 034502.
  * (158) Bonati, C.; D’Elia, M.; Mariti, M.; Martinelli, G.; Mesiti, M.; Negro, F.; Sanfilippo, F.; Villadoro, G. Axion phenomenology and $\theta$-dependence from $N_{f}=2+1$ lattice QCD. _JHEP_ 2016, _1603_ , 155.
  * (159) Vonk, T.; Guo, F.K.; Meißner, U.G. Aspects of the QCD $\theta$-vacuum. _JHEP_ 2019, _1906_ , 106.
  * (160) Vicari, E.; Panagopoulos, H. Theta dependence of SU(N) gauge theories in the presence of a topological term. _Phys. Rept._ 2009, _470_ , 93.
  * (161) Borsanyi, S.; Fodor, Z.; Guenther, J.; Kampert, K.-H.; Katz, S.D.; Kawanai, T.; Kovacs, T.G.; Mages, S.W.; Pasztor, A.; Pittler, F.; et al. Calculation of the axion mass based on high-temperature lattice quantum chromodynamics. _Nature_ 2016,_539_ , 69.


# HORTENSIA, a program package for the simulation of nonadiabatic
autoionization dynamics in molecules

Kevin Issler  Julius-Maximilians-Universität Würzburg, Institut für
Physikalische und Theoretische Chemie, Emil-Fischer-Str. 42, 97074 Würzburg,
Germany  Roland Mitrić<EMAIL_ADDRESS>Julius-Maximilians-
Universität Würzburg, Institut für Physikalische und Theoretische Chemie,
Emil-Fischer-Str. 42, 97074 Würzburg, Germany  Jens Petersen
<EMAIL_ADDRESS>Julius-Maximilians-Universität Würzburg,
Institut für Physikalische und Theoretische Chemie, Emil-Fischer-Str. 42,
97074 Würzburg, Germany

###### Abstract

We present a program package for the simulation of ultrafast vibration-induced
autoionization dynamics in molecular anions in the manifold of the adiabatic
anionic states and the discretized ionization continuum. This program, called
HORTENSIA ($\underline{Ho}$pping $\underline{r}$eal-time
$\underline{t}$rajectories for $\underline{e}$lectron-ejection by
$\underline{n}$onadiabatic $\underline{s}$elf-$\underline{i}$onization in
$\underline{a}$nions), is based on the nonadiabatic surface-hopping
methodology, wherein nuclei are propagated as an ensemble along classical
trajectories in the quantum-mechanical potential created by the electronic
density of the molecular system. The electronic Schrödinger equation is
numerically integrated along the trajectory, providing the time evolution of
electronic state coefficients, from which switching probabilities into
discrete electronic states are determined. In the case of a discretized
continuum state, this hopping event is interpreted as the ejection on an
electron. The derived diabatic and nonadiabatic couplings in the time-
dependent electronic Schrödinger equation are calculated from anionic and
neutral wavefunctions obtained from quantum chemical calculations with
commercially available program packages interfaced with our program.

Based on this methodology, we demonstrate the simulation of autoionization
electron kinetic energy spectra that are both time- and angle-resolved. In
addition, the program yields data that can be interpreted easily with respect
to geometric characteristics such as bonding distances and angles, which
facilitates the detection of molecular configurations important for the
autoionization process.

Moreover, useful extensions are included, namely generation tools for initial
conditions and input files as well as for the evaluation of output files both
through console commands and a graphical user interface.

††preprint: AIP/123-QED

##  I Introduction

After generation of a temporary molecular anion through electron attachment,
there are three possible competing relaxation mechanisms.Illenberger and
Momigny (1992) These are a) radiative deactivation, assuming that there is a
lower-lying anion state that is stable with respect to ionization, b)
dissociative electron attachment, in which the captured electron induces
geometric change in the molecule resulting in fragmentation into more stable
products, a neutral and an anionic subsystem. And lastly, c) autoionization,
in which after a finite period of time the metastable state decays via
electron ejection. The process of dissociative electron attachment is observed
for example in DNA, where capture of low-energy electrons leads to single and
double strand breaksBoudaïffa  _et al._ (2000); Martin  _et al._ (2004), or in
a variety of substances in nanoscale thin filmsArumainayagam  _et al._ (2010).
Prominent examples for autoionization include excited dipole- and quadrupole-
bound anions with binding energies slightly below the ionization thresholdLiu
_et al._ (2013, 2014); Zhu, Liu, and Wang (2017); Yuan, Zhang, and Wang
(2022), intermolecular Coulombic decay at the FADH- cofactor involved in DNA-
photolesion repairHarbach  _et al._ (2013) and autoionization induced by
vibrational excitation in organic moleculesAdams, Schneider, and Weber (2010);
Bull, West, and Verlet (2016); DeVine  _et al._ (2017, 2018); Adams, Hansen,
and Weber (2019); Anstöter  _et al._ (2020). Generally the finite lifetime of
a metastable state with respect to autoionization can vary strongly from only
a few femtosecondsNaff, Cooper, and Compton (1968); Jordan and Burrow (1987)
up to millisecondsNaff, Cooper, and Compton (1968); Suess, Parthasarathy, and
Dunning (2002). Recently, several experiments have provided insights into the
dynamics of such processes in dipole- and quadrupole-bound organic anions on a
(sub-)picosecond timescale.Bull, West, and Verlet (2016); Bull, Anstöter, and
Verlet (2019); Anstöter  _et al._ (2020); Verlet  _et al._ (2020); Kang, An,
and Kim (2020); Kang  _et al._ (2021)

Although the process of autoionization is well-known and -observed
experimentally by a multitude of methods, as can be seen in the references
given above, the theoretical description of autoionizing systems is
challengingHerbert (2015), especially if one is interested in the mechanistic
details of the intricate ultrafast relaxation dynamics. Autoionization
processes can follow different general mechanisms, depending on how energy is
redistributed among the system’s degrees of freedom. Besides a purely
electronic variant, where already the electronic energy of the system lies
above the ionization threshold and electron ejection may proceed via
tunneling, there is also the possibility of a nonadiabatic mechanism in which
rotational or vibrational energy of the nuclei is transformed into the kinetic
energy of the ejected electron.

In the following, we focus on the case of vibrational autoionization. This
process can thus be viewed as a nonadiabatic transition between a
vibrationally excited bound N-electron system and continuum electronic states
consisting of an N-1 electron molecular core and a free electron. Early
theoretical treatments have focused on the computation of ionization
ratesBerry (1966); Acharya, Kendall, and Simons (1984); Simons (1999) as well
as on establishing propensity rules for the ionization transitionsSimons
(1981). While a full dynamical treatment of vibrational autoionization is
highly desirable, an entirely quantum-dynamical approach is computationally
prohibitive. As an alternative, a mixed quantum-classical ansatz can be
considered, further motivated by the success of this type of methodology in
the description of bound-state nonadiabatic processes and the simulation of
time-resolved spectroscopic signals.Humeniuk  _et al._ (2013); Persico and
Granucci (2014); Richings  _et al._ (2015); Crespo-Otero and Barbatti (2018);
Röder  _et al._ (2019) Although to date there have been several
implementations of mixed quantum-classical dynamics simulations for bound-
state problems made publicly availableBarbatti  _et al._ (2014); Mai,
Marquetand, and González (2018); Du and Lan (2015), no program addressing the
simulation of vibration-induced autoionization processes has been published so
far.

Therefore, in this work we present the program package implementing our
approach to describe vibrational autoionization through quantum-classical
dynamics in the framework of the surface-hopping methodology in the manifold
of bound and continuum electronic states as described recentlyIssler, Mitrić,
and Petersen (2023a). Therein, nuclear motion is considered classically, while
the electronic system is treated quantum-mechanically. Nonadiabatic
transitions between electronic states accompanied by change of the classical
vibrational energy of the molecule describe the energy exchange between the
two subsystems. With this program package and the underlying methodology, one
is able to gain insight into the geometric and electronic evolution in the
course of the autoionization process as well as to calculate the time-,
energy- and angle-distribution of the generated free electrons, which serve as
experimental observables for monitoring autoionization dynamics.

We illustrate our program on the example of the 2-cyanopyrrolide anion, which
bears a dipole-bound excited state slightly below the electron detachment
threshold while the vibrationally excited states are metastable and decay via
autoionization.Yuan, Zhang, and Wang (2022)

In the following section a brief theoretical description of the method is
given. In section III an overview of the actual implementation is provided.
The subsequent section IV details performance-related issues, namely quality
of approximations in the theory and runtime and memory optimization within the
program, as well as a dynamics simulation example for the 2-cyanopyrrolide
anion. Finally in section V a conclusion and outlook are given.

##  II Theory

Our methodological framework is based on the surface-hopping procedure as
proposed by TullyTully (1990), in which the coupled electron-nuclear dynamics
of molecular systems is approached in a quantum-classical fashion.
Specifically, the nuclei are propagated classically according to Newton’s
equations of motion,

$\displaystyle
M\ddot{\textbf{R}}(t)=\mathbf{F}_{i}(\mathbf{R}[t])\equiv-\nabla_{\textbf{R}}E_{i}(\textbf{R}[t]),$
(1)

where the force $\mathbf{F}_{i}(\mathbf{R}[t])$ is obtained as the negative
gradient of the electronic potential energy surface (PES)
$E_{i}(\textbf{R}[t])$. In the above equation, $M$ denotes a diagonal matrix
containing the nuclear masses. For an ensemble of initial conditions, this
leads to trajectories $\textbf{R}(t)$ moving on the given PES. Simultaneously,
the electronic time-dependent Schrödinger equation

$\displaystyle
i\hbar\dot{\Psi}(\textbf{r};\textbf{R}[t])=\hat{H}_{el}\Psi(\textbf{r};\textbf{R}[t]),$
(2)

with the electronic Hamiltonian $\hat{H}_{el}$ is solved. The electronic
wavefunction can be expanded into a set of orthonormal basis states, which in
the case of autoionization includes bound states $\Phi_{m^{\prime}}$ (denoted
with a primed index) as well as continuum states
$\tilde{\Phi}_{m^{\prime\prime}}$ (denoted with a double-primed index):

$\displaystyle\Psi\big{(}\textbf{r},\textbf{R}[t],t\big{)}=$
$\displaystyle\sum_{m^{\prime}}c_{m^{\prime}}(t)\Phi_{m^{\prime}}\big{(}\textbf{r},\textbf{R}[t]\big{)}\
+$ $\displaystyle\sum_{m^{\prime\prime}}\int\\!d^{3}\textbf{k}\
\tilde{c}_{m^{\prime\prime}}(\textbf{k},t)\tilde{\Phi}_{m^{\prime\prime}}(\textbf{k},\textbf{r},\textbf{R}[t]),$
(3)

where k denotes the continuously varying wave vector of the free electron,
while $m^{\prime\prime}$ is the quantum number of the remaining neutral state.
We assume the wavefunctions $\Phi_{m^{\prime}}$ and
$\tilde{\Phi}_{m^{\prime\prime}}$ to be single Slater determinants (ground
state) or an expansion of singly excited Slater determinants (excited state).
In the frame of the presented methodology we discretize the continuum states,
leading to

$\displaystyle\int
d^{3}\textbf{k}\,\tilde{c}_{m^{\prime\prime}}(\textbf{k},t)\tilde{\Phi}_{m^{\prime\prime}}(\textbf{k},\textbf{r},\textbf{R}[t])$
$\displaystyle\approx\sum_{i}(\Delta{\cal
V}_{k})^{\frac{1}{2}}\tilde{c}_{m^{\prime\prime}}(\textbf{k}_{i},t)(\Delta{\cal
V}_{k})^{\frac{1}{2}}\tilde{\Phi}_{m^{\prime\prime}}(\textbf{k}_{i},\textbf{r},\textbf{R}[t])$
$\displaystyle\approx\sum_{i}c_{m^{\prime\prime}}(\textbf{k}_{i},t)\Phi_{m^{\prime\prime}}(\textbf{k}_{i},\textbf{r},\textbf{R}[t]),$
(4)

where $\Delta{\cal V}_{k}$ denotes the volume element in $k$-space and the
discretized and continuum state expansion coefficients are related according
to $c_{m^{\prime\prime}}(\textbf{k}_{i},t)=(\Delta{\cal
V}_{k})^{\frac{1}{2}}\tilde{c}_{m^{\prime\prime}}(\textbf{k}_{i},t)$. The
actual determination of the wave vectors and the implementation of the
discretization procedure are explained in detail in the next chapter.

Insertion of Eq. (3) into the time-dependent Schrödinger equation (2),
multiplication from the left by an eigenstate $\bra{\Phi_{n}}$ and evaluation
of the arising terms leads to a set of coupled differential equations for the
electronic state coefficients $c_{n}$:

$\displaystyle\dot{c}_{n}(t)=\sum_{j}\left[-\frac{i}{\hbar}H_{nm}(\textbf{R}[t])-D_{nm}(\textbf{R}[t])\right]c_{m}(t),$
(5)

with the matrix elements of the electronic Hamiltonian
$H_{nm}=\braket{\Phi_{n}}{H_{el}}{\Phi_{m}}$ and the nonadiabatic couplings
$D_{nm}=\braket{\Phi_{n}}{\dot{\Phi}_{m}}=\dot{\textbf{R}}\cdot\braket{\Phi_{n}}{\nabla_{R}}{\Phi_{m}}$.
These can be divided into separate expressions for the bound and continuum
states, resulting in the diabatic and nonadiabatic couplings between two bound
anion states,

$\displaystyle H_{n^{\prime}m^{\prime}}$
$\displaystyle=\braket{\Phi_{n^{\prime}}}{\hat{H}}{\Phi_{m^{\prime}}}$ (6)
$\displaystyle D_{n^{\prime}m^{\prime}}$
$\displaystyle=\braket{\Phi_{n^{\prime}}}{\dot{\Phi}_{m^{\prime}}},$ (7)

and between a bound and a discretized continuum state,

$\displaystyle H_{n^{\prime\prime}m^{\prime}}(\textbf{k}_{i})$
$\displaystyle=(\Delta{\cal
V}_{k})^{\frac{1}{2}}\braket{\tilde{\Phi}_{n^{\prime\prime}}(\textbf{k}_{i})}{\hat{H}}{\Phi_{m^{\prime}}}$
(8) $\displaystyle D_{n^{\prime\prime}m^{\prime}}(\textbf{k}_{i})$
$\displaystyle=\braket{\Phi_{n^{\prime\prime}}(\textbf{k}_{i})}{\dot{\Phi}_{m^{\prime}}}=(\Delta{\cal
V}_{k})^{\frac{1}{2}}\braket{\tilde{\Phi}_{n^{\prime\prime}}(\textbf{k}_{i})}{\dot{\Phi}_{m^{\prime}}}.$
(9)

In the above equations, the approximation to neglect the coupling terms
between the continuum states has been introduced. The discretized continuum
states consist of an antisymmetrized product of a bound $N-1$ electron neutral
state and a molecular scattering state of the free electron

$\displaystyle\tilde{\Phi}_{n^{\prime\prime}}(\textbf{k}_{i})={\cal
A}\left(\Phi^{\textrm{(n)}}_{n^{\prime\prime}}\cdot\psi(\textbf{k}_{i})\right).$
(10)

The simplest approximation to the free electron states in the presence of a
neutral molecular core are plane waves

$\displaystyle\psi(\textbf{k}_{i})\approx{\cal
N}\textrm{e}^{i\textbf{k}_{i}\cdot\textbf{r}}$ (11)

with a normalization constant ${\cal N}=(2\pi)^{-3/2}$ to satisfy the
orthonormality demanded in Eq. (3). Since this function would be completely
independent on the electronic and nuclear configuration of the molecular core,
which is a strong simplification, the plane waves are orthogonalized with
respect to the anion’s molecular orbitals (MOs) $\phi_{m}$ to include (at
least to a certain degree) dependence on the molecular structure according to

$\displaystyle\tilde{\psi}(\textbf{k}_{i})$ $\displaystyle=(2\pi)^{-3/2}{\cal
N}_{ortho}\left(\textrm{e}^{i\textbf{k}_{i}\cdot\textbf{r}}-\sum_{m}^{\mathrm{occ}}\braket{\phi_{m}}{\textrm{e}^{i\textbf{k}_{i}\cdot\textbf{r}}}\phi_{m}\right)$
$\displaystyle={\cal
N}_{ortho}\left(\psi(\textbf{k}_{i})-\sum_{m}^{\mathrm{occ}}\braket{\phi_{m}}{\psi(\textbf{k}_{i})}\,\phi_{m}\right),$
(12)

with the normalization constant

$\displaystyle{\cal
N}_{ortho}=\left(1-\sum_{m}^{\mathrm{occ}}\big{|}\braket{\phi_{m}}{\psi(\textbf{k}_{i})}\big{|}^{2}\right)^{-1/2}$
(13)

arising from the orthogonalization.

Notably, the summation over $m$ includes the occupied MOs in all ’relevant’
Slater determinants of all considered electronic states, that is, we
considered all determinants which are needed to sufficiently represent the
ground state and full CIS wavefunction of the excited state. Beginning from
the highest contribution to a wavefunction, determinants are included until a
specific percentage or a user-adjusted maximum number of configurations per
electronic state is reached (95 % / 10 configurations in the case of
vinylideneIssler, Mitrić, and Petersen (2023a)). Considering for now the
special case where only the anion’s ground state is included, the used MOs are
simply the energetically lowest ones up to the highest-occupied molecular
orbital (HOMO).

The overlap integral between a plane wave and an MO present in Eq. (13),
$\braket{\phi_{m}}{\psi(\textbf{k}_{i})}$, can be computed analytically by
expanding the MO into the Gaussian atomic orbital (AO) basis, with the
integral involving a single AO $|\nu\rangle$ given by

$\displaystyle\braket{\nu}{\psi(\textbf{k})}=$ $\displaystyle\
(2\pi)^{-3/2}\int
d^{3}\mathbf{r}\,\textrm{e}^{i\textbf{k}\cdot\textbf{r}}\varphi_{\nu}(\textbf{r})$
$\displaystyle=$ $\displaystyle\
(2\alpha_{\nu})^{-3/2}\exp{\left(i\textbf{k}\cdot\textbf{A}_{\nu}-\frac{k^{2}}{4\alpha_{\nu}}\right)}$
$\displaystyle\times\prod_{j=x,y,z}(-i\sqrt{4\alpha_{\nu}})^{-n_{\nu,j}}H_{n_{\nu,j}}\left(\frac{k_{j}}{\sqrt{4\alpha_{\nu}}}\right),$
(14)

where the $H_{n_{\nu,j}}$ are the Hermite polynomials of order $n_{\nu,j}$.

###  II.1 Electronic coupling terms

There are anionic systems, for example the vinylidene anionIssler, Mitrić, and
Petersen (2023a), that do not support a bound excited state, in which case the
consideration of only the ground state and the continuum in the process of
autoionization is sufficient. Besides that, for example in molecules
exhibiting dipole-bound excited states Jordan and Wang (2003); Yuan, Zhang,
and Wang (2022); Issler, Mitrić, and Petersen (2023b), several bound anionic
states and the interaction among them are relevant as well. Nonetheless, to
keep the formalism concise, if not noted otherwise we discuss in the following
the electronic coupling terms for the special case of both anion and neutral
molecule being in their respective electronic ground states, which in turn are
represented by a single Slater determinant. The generalization to excited
states and/or multideterminantal wavefunctions is straightforward.Issler,
Mitrić, and Petersen (2023b) We denote the bound anionic ground state
wavefunction by $\ket{\Phi_{0}}$ and the continuum wavefunctions by
$\ket{\Phi_{i}}$, the latter being constructed as an antisymmetrized product
of the neutral ground state and a free electron state function with wave
vector $\textbf{k}_{i}$, similar to Eq. (10).

####  II.1.1 Diabatic couplings

In the case of two adiabatic bound anion states, the coupling matrix elements
$H_{n^{\prime}m^{\prime}}$ given in Eq. (6) yield zero for all $n^{\prime}\neq
m^{\prime}$ since these states are orthonormal eigenstates of the electronic
Hamiltonian.

On the other hand, since in our methodology the bound and continuum state
wavefunctions are constructed using separate quantum-chemical calculations for
the anion and neutral, and the free electron wavefunction is taken as a plane
wave, the continuum state functions are crude approximations to the actual
adiabatic eigenfunctions of the electronic Hamiltonian for the $N$-electron
system and therefore, diabatic couplings between the bound and continuum
electronic states arise.

As elaborated in detail in Ref. Issler, Mitrić, and Petersen, 2023a, according
to Eq. (8) and defining $V_{i0}^{\mathrm{dia}}(\textbf{k}_{i})$ as

$\displaystyle
H_{i0}(\textbf{k}_{i})\equiv\braket{\Phi_{i}}{\hat{H}}{\Phi_{0}}\equiv(\Delta{\cal
V}_{k})^{\frac{1}{2}}\,V^{\mathrm{dia}}_{i0}(\textbf{k}_{i}),$ (15)

the diabatic coupling between a bound and a continuum state can be written in
terms of the AO basis as

$\displaystyle V^{\mathrm{dia}}_{i0}(\textbf{k}_{i})$ $\displaystyle=$
$\displaystyle\sum_{\lambda\mu\nu}\Bigg{[}$ $\displaystyle
A_{\lambda\mu\nu}\Big{(}\braket{\mathbf{k}_{i}\lambda}{\mu\nu}-\sum_{\sigma}B_{\sigma}\braket{\sigma\lambda}{\mu\nu}\Big{)}+$
$\displaystyle\bar{A}_{\lambda\mu\nu}\Big{(}\braket{\mathbf{k}_{i}\lambda}{\mu\nu}-\sum_{\sigma}B_{\sigma}\braket{\sigma\lambda}{\mu\nu}\Big{)}\Bigg{]}.$
(16)

In this formula the Greek letters denote the AO basis functions,
$\braket{\mathbf{k}_{i}\lambda}{\mu\nu}$ is an electron-electron repulsion
integral and
$\braket{\mathbf{k}_{i}\lambda}{\mu\nu}=\braket{\textbf{k}_{i}\lambda}{\mu\nu}-\braket{\textbf{k}_{i}\lambda}{\nu\mu}$
its antisymmetrized variant. The prefactors $A_{\lambda\mu\nu}$,
$\bar{A}_{\lambda\mu\nu}$ and $B_{\sigma}$ comprise AO expansion coefficients
and overlap integrals and are defined as follows (assuming that the extra
electron of the anion has $\alpha$ spin):

$\displaystyle A_{\lambda\mu\nu}=$
$\displaystyle\sum_{n}^{\mathrm{occ},\alpha}\sum_{q,p<q}^{\mathrm{occ},\alpha}(-1)^{n+p+q-1}\mathrm{det}\
\mathbf{S}_{in,pq}$
$\displaystyle\times\left(c_{\lambda}^{(n)}-\sum_{u}^{\mathrm{occ},\alpha}c_{\lambda}^{(u)}S_{nu}\right)c_{\mu}^{(p)}c_{\nu}^{(q)}$
(17) $\displaystyle\bar{A}_{\lambda\mu\nu}=$
$\displaystyle\sum_{\bar{n}}^{\mathrm{occ},\beta}\sum_{p}^{\mathrm{occ},\alpha}\sum_{\bar{q}}^{\mathrm{occ},\beta}(-1)^{\bar{n}+p+\bar{q}-1}\mathrm{det}\
\mathbf{S}_{i\bar{n},p\bar{q}}$
$\displaystyle\times\left(c_{\lambda}^{(\bar{n})}-\sum_{\bar{u}}^{\mathrm{occ},\beta}c_{\lambda}^{(\bar{u})}S_{\bar{n}\bar{u}}\right)c_{\mu}^{(p)}c_{\nu}^{(\bar{q})}$
(18) $\displaystyle B_{\sigma}=$
$\displaystyle\sum_{r}^{\mathrm{occ},\alpha}\sum_{\rho}c_{\sigma}^{(r)}c_{\rho}^{(r)}\braket{\textbf{k}_{i}}{\rho},$
(19)

where the indices (including their variants with an overbar) $p,q,r$ refer to
anion MOs, $n,u$ to neutral MOs, and $\mathrm{det}\ \mathbf{S}_{in,pq}$
denotes the minor determinant of the overlap matrix between continuum and
bound state orbitals where the rows of the free electron orbital
$\tilde{\psi}(\mathbf{k}_{i})$ and neutral orbital $\chi_{n}$ as well as the
columns of anion orbitals $\phi_{p}$ and $\phi_{q}$ have been deleted. For the
full derivation of these equations the reader is referred to Ref. Issler,
Mitrić, and Petersen, 2023a.

####  II.1.2 Nonadiabatic couplings

The nonadiabatic coupling terms as defined in Eqs. (7) and (9) are calculated
using the finite-difference approximation for the time derivative, which leads
to

$\displaystyle D_{i0}(t)$
$\displaystyle=\braket{\Phi_{i}(t)}{\frac{d}{dt}{\Phi}_{0}(t)}$ (20)
$\displaystyle\approx\frac{1}{2\Delta t}\Big{(}\braket{\Phi_{i}(t-\Delta
t)}{\Phi_{0}(t)}-\braket{\Phi_{i}(t)}{\Phi_{0}(t-\Delta t)}\Big{)}$ (21)

In the case of two anionic bound states, these terms are evaluated according
to Refs. Mitrić, Werner, and Bonačić-Koutecký, 2008; Werner  _et al._ , 2008;
Werner, Mitrić, and Bonačić-Koutecký, 2010.

One can simplify the arising terms by integrating over all but one electron
coordinate. For the first term of Eq. (21) this yields

$\displaystyle\braket{\Phi_{i}(t^{\prime})}{\Phi_{0}(t)}=N^{-1/2}\braket{\tilde{\psi}(\textbf{k}_{i},t^{\prime})}{\psi^{D}(t^{\prime},t)},$
(22)

where we have abbreviated $t^{\prime}=t-\Delta t$ and have defined the one-
electron function $\psi^{D}(t^{\prime},t)$, which is an analog to a molecular
Dyson orbital with the $N$\- and $N-1$\- wavefunctions taken at different time
steps and geometries. Using Eqs. (12) and (22) the resulting nonadiabatic
coupling terms read

$\displaystyle D_{i0}(\textbf{k}_{i},t)=\frac{(\Delta{\cal
V}_{k})^{\frac{1}{2}}{\cal N}_{ortho}}{2\sqrt{N}\Delta
t}\Big{[}\braket{\psi(\textbf{k}_{i})}{\psi^{D}(t^{\prime},t)}-\braket{\psi(\textbf{k}_{i})}{\psi^{D}(t,t^{\prime})}$
$\displaystyle-\sum_{n}\braket{\psi(\textbf{k}_{i})}{\phi_{n}(t)}\,\braket{\phi_{n}(t^{\prime})}{\psi^{D}(t^{\prime},t)}$
$\displaystyle+\sum_{n}\braket{\psi(\textbf{k}_{i})}{\phi_{n}(t)}\,\braket{\phi_{n}(t)}{\psi^{D}(t,t^{\prime})}\Big{]}.$
(23)

###  II.2 Adiabatic ionization and electronic decay

The main focus of the above presented methodology lies on describing the
nonadiabatic process of vibrational autoionization. However, in the course of
the molecule’s dynamical evolution instances can occur where the occupied
anionic state becomes unbound as the result of changes in nuclear geometry. In
this case, ionization is possible as an exclusively adiabatic electronic
process without coupling to the nuclear motion. This process can be included
approximately in our method by simulating the temporal spread of the ejected
electron as a wavepacket evolving freely in space. As a quantitative measure,
the electronic spatial extent, i.e., the expectation value of
$\hat{\mathbf{r}}^{2}$, is calculated as a function of time.

Specifically, once a time step is reached where the VDE has become negative,
the highest-occupied orbital of the last bound geometry,
$\phi(\textbf{r},t_{0})$, is used as the initial free electronic wavepacket.
In the case where one only considers the anionic ground state, this
corresponds to the HOMO. If also an excited state is involved, natural
transition orbitals (NTOs)Martin (2003) are calculated and the highest-
occupied and lowest-unoccupied NTO (HONTO and LUNTO) are used for the anionic
ground and excited state, respectively. Such an electronic wavepacket is then
propagated in time and its spatial extent is evaluated according to

$\displaystyle\braket{\hat{\mathbf{r}}^{2}}(t)$
$\displaystyle=\braket{\phi(\mathbf{r},t)}{\hat{\mathbf{r}}^{2}}{\phi(\mathbf{r},t)}$
$\displaystyle=\sum_{\mu\nu}c_{\mu}c_{\nu}\braket{\varphi_{\mu}(\mathbf{r},t)}{\hat{\mathbf{r}}^{2}}{\varphi_{\nu}(\mathbf{r},t)}.$
(24)

Here $\varphi_{\mu,\nu}$ denote the Gaussian atomic basis functions freely
propagated in time:

$\varphi_{\mu}(\mathbf{r},t)=\int
d^{3}\mathbf{r}^{\prime}\,K(\mathbf{r},\mathbf{r}^{\prime},t,0)\varphi_{\mu}(\mathbf{r}^{\prime},0)$
(25)

with the free electron propagator

$K(\mathbf{r},\mathbf{r}^{\prime},t,0)=\Braket{\mathbf{r}}{\mathrm{e}^{-i\hat{\mathbf{p}}^{2}t/2m_{e}\hbar}}{\mathbf{r}^{\prime}}.$
(26)

Using Cartesian Gaussian basis functions of $s$, $p$ and $d$ type one obtains
the following analytic expression for the electronic wavepacket:

$\displaystyle\varphi_{\mu}(\mathbf{r},t)=N_{l_{x}l_{y}l_{z}}\mathrm{e}^{-\frac{\alpha}{1+i\beta
t}\mathrm{r}^{2}}\left[-\Lambda\frac{i\beta t}{2\alpha}(1+i\beta
t)^{-\frac{5}{2}}+\right.$ $\displaystyle\left.(1+i\beta
t)^{-\frac{3}{2}-\sum_{j}l_{j}}\prod_{j=x,y,z}(r_{j}-A_{j})^{l_{j}}\right],$
(27)

where A is the spatial center of the respective basis function, $l_{i}$
denotes the angular momentum quantum number for the $i$’th spatial direction
and $\Lambda$ is a constant that is unity if one of the $l_{i}=2$ and zero if
all $l_{i}<2$. The AO integrals in Eq. (24) are calculated with an
implementation of the McMurchie-Davidson schemeMcMurchie and Davidson (1978).
To relate the spatial extent in a simple way to the lifetime of the unbound
state, an auxiliary spherically symmetric electron distribution is considered
which within the initially determined radius
$r_{0}=\sqrt{\braket{\textbf{r}^{2}}(t_{0})}$ contains a probability of 99%.
Subsequently, with $\braket{\textbf{r}^{2}}$ increasing with time, the
probability within $r_{0}$ decreases, giving rise to a population decay curve
which can be related to a time constant $\tau$. The latter is incorporated
into the propagation of the electronic wavefunction given by Eq. (5) by adding
an imaginary component to the electronic state energy,

$E^{\mathrm{(a)}}\rightarrow E^{\mathrm{(a)}}-\frac{i\hbar}{2\tau},$ (28)

which leads to an exponential population decay due to adiabatic ionization in
regions where the VDE is negative for the given electronic state.

###  II.3 Surface-hopping procedure

Solution of the set of Eqs. (5) along a nuclear trajectory yields the time-
dependent electronic state coefficients $c_{n}(t)$. Within the surface-hopping
methodology, a switch from the occupied bound electronic state $n$ to any
other state $m$ is determined by the hopping probability depending on the
electronic state populations $\rho_{nn}=|c_{n}|^{2}$, which is

$\displaystyle P_{n\rightarrow
m}=-\frac{\dot{\rho}_{nn}}{\rho_{nn}}\frac{\dot{\rho}_{mm}}{\sum_{k}\dot{\rho}_{kk}}\Delta
t$ (29)

for $\dot{\rho}_{nn}<0$ and $\dot{\rho}_{mm}>0$ and zero in any other
instance. In the above expression, the sum over $k$ includes all states with
$\dot{\rho}_{kk}>0$. In case a surface hop occurs, to ensure energy
conservation the nuclear velocities are rescaled such that for kinetic
energies $T$ and electronic potential energies $E_{n}$ of anion (a) and
neutral (n) the following conditions are fulfilled:

$\displaystyle
T^{\prime\textrm{(a)}}=T^{\textrm{(a)}}+E_{n}^{\textrm{(a)}}-E_{m}^{\textrm{(a)}}$
(30)

for a hop between anionic bound states and

$\displaystyle
T^{\prime\textrm{(n)}}=E_{n}^{\textrm{(a)}}+T^{\textrm{(a)}}-E_{m}^{\textrm{(n)}}-E_{\textrm{el}}(\textbf{k}_{i})$
(31)

for a hop into the continuum (i.e. autoionization). For a more detailed
description of the hopping procedure the reader is referred to Ref. Mitrić,
Petersen, and Bonačić-Koutecký, 2011.

##  III Program implementation

Figure 1: Schematic of the dynamics procedure as implemented in the HORTENSIA
program package. The box coloration matches the specific tasks to the program
modules as follows: red: wigner/wignerEnsemble.py, blue: external quantum-
chemistry program, orange: nuclearDynamics.py, green: populationDynamics.py

In the following chapter a detailed account of how the theory is actually
implemented in the program package will be provided. For an easier
understanding, in Fig. 1 the program flow is displayed schematically, with a
color code indicating the module handling the respective task.

Starting from the generation of an ensemble of nuclear coordinates
$\textbf{R}(t)$ and velocities $\dot{\textbf{R}}(t)$ at the time
$t=t_{initial}$ using the wignerEnsemble module in the wigner folder (red), a
first quantum-chemical calculation is performed by an external quantum-
chemistry program - to date these include Gaussian09/Gaussian16 Frisch  _et
al._ (2016) and QChem Shao  _et al._ (2015) (blue) - which yields the forces
from which the accelerations $\ddot{\textbf{R}}(t)$ of the nuclei are
computed. The nuclei are then propagated by integration of Newton’s equations
of motion for one nuclear time step using the nuclearDynamics module (orange).
With the new nuclear coordinates $\textbf{R}(t+\Delta t)$, a new set of
quantum-chemical calculations can be performed, yielding the new energy
gradients necessary for the evaluation of the velocities
$\dot{\textbf{R}}(t+\Delta t)$. With the quantum-chemical calculations at $t$
and $t+\Delta t$, one is now able to construct the electronic continuum states
as well as the coupling matrices of the diabatic and nonadiabatic couplings
using the populationDynamics module (green). From this point, the electronic
state coefficients $\textbf{c}(t)$ are propagated in parallel to the nuclear
dynamics by integrating the electronic Schrödinger equation, yielding
$\textbf{c}(t+\Delta t)$. These are utilized to compute hopping probabilities
from the occupied bound state to all other (bound and continuum) states. The
switching between the states is induced stochastically according to the
respective hopping probabilities given in Eq. (29). After writing the results
into the various output files time is shifted to $t=t+\Delta t$, thereby
completing one time step.

To make this initial overview more specific, in the following the underlying
algorithms are explained in more detail.

###  III.1 Electronic structure calculation

All electronic structure and energy gradient calculations can be performed by
using any Kohn-Sham (TD)-DFT level of theory provided within the Gaussian09,
Gaussian16 or QChem program packages. The AO basis set needs to be defined
explicitly in a separate input file, thus also allowing for additional
augmentation of basis sets, which is of utmost importance when describing
molecular anions.Herbert (2015); Skurski, Gutowski, and Simons (2000) The
handlerG09 and handlerQChem modules provide an interface to the external
programs by creating input files and calling the respective programs. The
dysonG09 and dysonQChem modules contain classes that parse the external output
files and organize the data into the form needed in the program.

###  III.2 Generation of initial conditions

The initial nuclear coordinates and velocities are determined by stochastic
sampling of an appropriate probability distribution function for the harmonic
normal modes of the system. These can be computed from the electronic Hessian
matrix at an optimized geometry of the studied molecule. For molecules in the
vibrational ground state as well as for a thermal ensemble of molecules, the
Wigner function

$\displaystyle\rho_{W}(\\{Q_{i},P_{i}\\})=\frac{1}{(\pi\hbar)^{N}}\prod_{i=1}^{N}\alpha_{i}(T)\,\mathrm{exp}\left(-\frac{\alpha_{i}(T)}{\hbar\omega_{i}}(P_{i}^{2}+\omega_{i}^{2}Q_{i}^{2})\right)$
(32)

with

$\displaystyle\alpha_{i}(T)=\tanh{\bigg{(}\frac{\hbar\omega_{i}}{2k_{B}T}\bigg{)}}$
(33)

is employed, where $\\{Q_{i},P_{i}\\}$ denote the normal coordinates and
momenta, $\omega_{i}$ is the angular frequency of normal mode $\nu_{i}$ and
$T$ the thermodynamic temperature.

Besides these cases, in experiments investigating vibration-induced
autoionization another type of initial conditions is often important in which
one or more normal vibrations of the system are excited by laser irradiation.
In principle, the respective initial conditions could be also generated by
using a Wigner function. However, Wigner functions for excited vibrational
states can assume negative values and can thus not be directly identified with
a probability distribution. A possible approach might be to regard the
positive and negative parts of the Wigner function separately as probability
distributions and to run a "positive" and a "negative" ensemble of initial
conditions, the final properties of the system then being obtained by
appropriate averaging. As a more efficient alternative, which gets on with
only one single ensemble, we employ a positive definite probability
distribution constructed from the excited-vibrational state wavefunctions in
position and momentum space,

$\displaystyle\rho^{(i)}_{\upsilon}(Q_{i},P_{i})=|\chi^{(i)}_{\upsilon}(Q_{i})|^{2}|\tilde{\chi}^{(i)}_{\upsilon}(P_{i})|^{2},$
(34)

where $\chi^{(i)}_{\upsilon}(Q_{i})$ and
$\tilde{\chi}^{(i)}_{\upsilon}(P_{i})$ are the harmonic oscillator
wavefunctions for quantum state $\upsilon$ of normal mode $\nu_{i}$ in
position and momentum space, respectively.

###  III.3 Nuclear dynamics

Given Newton’s equations of motion (1), the nuclei are propagated by numerical
solution using the velocity Verlet algorithm Swope  _et al._ (1982) for a
user-defined time step. Within this algorithm, the nuclear coordinates at
$t+\Delta t$ are obtained from a Taylor series expansion around the
coordinates at $t$:

$\displaystyle\textbf{R}(t+\Delta
t)\approx\textbf{R}(t)+\dot{\textbf{R}}(t)\Delta
t+\frac{1}{2}M^{-1}\textbf{F}(t)\Delta t^{2},$ (35)

where in the last term the acceleration has been formulated using the force F
given by the electronic potential energy gradient (cf. Eq. (1)). With the new
nuclear coordinates, the force at $t+\Delta t$ can be evaluated, giving rise
to the new nuclear velocities

$\displaystyle\dot{\textbf{R}}(t+\Delta t)=\dot{\textbf{R}}(t)+\frac{\Delta
t}{2}M^{-1}\left[\textbf{F}(t)+\textbf{F}(t+\Delta t)\right].$ (36)

Due to the approximative nature of the algorithm above and the accuracy of the
calculated energy gradients, it is possible that the velocities develop small
overall translational or rotational components although the initial conditions
were determined with these degrees of freedom set at rest. These numerical
inaccuracies are detected, in the case of translational velocity by the shift
of the center of mass away from the origin of the coordinate system, in the
case of rotation by the calculation of the angular velocity according to

$\displaystyle\bm{\omega}_{rot}=I^{\,-1}\textbf{L}$ (37)

with the moment of inertia $I$ and the angular momentum L. The translational
and rotational portions of the nuclear velocities are then subtracted from the
total velocity and the remaining vibrational part is rescaled to ensure energy
conservation.

After each nuclear dynamics step, the new nuclear coordinates and velocities
are written into separate output files, the coordinates in a format of
consecutive xyz files which can be visualized easily by external software (for
example with the VMD program package Humphrey, Dalke, and Schulten (1996),
which is warmly recommended).

###  III.4 Electronic dynamics

Since the evaluation of electronic coupling terms in Eq. (5) is, apart from
the external quantum-chemistry calculations, the computationally most
expensive step in the dynamics, several approximations need to be implemented,
which will be discussed in the following

####  III.4.1 Calculation of coupling terms

Before calculating the coupling terms, the discretization procedure for the
generation of wave vectors needed to construct the continuum state
wavefunctions will be discussed. To uniformly discretize angular orientation
and kinetic energy of ejected electrons, it is natural to discretize angular
and energetic distribution separately. Since the kinetic energy of a plane
wave is

$\displaystyle
E_{kin}(\textbf{k}_{i})=\frac{\hbar^{2}|\textbf{k}_{i}|^{2}}{2m_{e}}$ (38)

and therefore proportional to the length of the wave vector squared, this
length is discretized such that the desired energy range is covered evenly.
For a given energy, the vector orientations are approximately evenly
distributed according to the Fibonacci sphere algorithm Swinbank and James
Purser (2006). The volume elements $\Delta{\cal V}_{k}$ needed for calculating
the bound-continuum couplings in Eqs. (8) and (9) are constructed as the
difference of spherical caps around the corresponding wave vectors with the
base diameter as an average over the six nearest points on the sphere
surrounding the vector.

In the diabatic coupling terms in the AO basis (Eq. (16)) two types of four-
center integrals are present: (i) such involving four Gaussian-type atomic
orbitals (GTOs), $\braket{\sigma\lambda}{\mu\nu}$. These are evaluated by
using the libcint library Sun (2015) within the PySCF program package Sun  _et
al._ (2018, 2020). (ii) integrals involving a plane wave of wave vector
$\mathbf{k}_{i}$ and three GTOs, $\braket{\mathbf{k}_{i}\lambda}{\mu\nu}$.
These terms can in principle be calculated analytically as outlined, e.g., in
Ref. Colle, Fortunelli, and Simonucci, 1987, but this is computationally
unfeasible for the present purpose since an immense number of plane waves has
to be included for a proper discretization of the ionization continuum.
Instead, the plane waves are approximated by their Taylor expansion around the
center of basis function $\ket{\mu}$, $\textbf{R}_{\mu}$. As will be discussed
in the Performance Section later on, for sufficient accuracy in the
approximation it is necessary to include not only the zero’th order term
(assuming the plane wave to be constant in the vicinity of the molecule), but
also the first-order term, resulting in the approximation

$\displaystyle\mathrm{e}^{i\textbf{k}\cdot\textbf{r}}$
$\displaystyle=\mathrm{e}^{i\textbf{k}\cdot\textbf{R}_{\mu}}\mathrm{e}^{i\textbf{k}\cdot(\textbf{r}-\textbf{R}_{\mu})}$
$\displaystyle\approx\mathrm{e}^{i\textbf{k}\cdot\textbf{R}_{\mu}}\left[1+i\textbf{k}\cdot(\textbf{r}-\textbf{R}_{\mu})\right].$
(39)

This leads to two terms for the two-electron integrals as follows:

$\displaystyle\braket{\mathbf{k}_{i}\lambda}{\mu\nu}\approx\mathrm{e}^{i\textbf{k}\cdot\textbf{R}_{\mu}}\left[\braket{\lambda}{\mu\nu}+i\textbf{k}\braket{\lambda}{\tilde{\mu}\nu}\right].$
(40)

In the above expression, $\ket{\tilde{\mu}}$ is an AO basis function with an
angular momentum quantum number by one higher than $\ket{\mu}$ while having
the same Gaussian exponent. This heavily reduces the amount of two-electron
integrals to be computed from $n_{AO}^{3}n_{PW}$ to
$n_{AO}^{2}[n_{AO}+n^{\prime}_{AO}]$, with $n_{AO}$ being the total number of
AO basis functions, $n^{\prime}_{AO}$ the total number of basis functions with
increased quantum number and $n_{PW}$ the total number of plane waves. For
instance, in the case of vinylidene in Ref. Issler, Mitrić, and Petersen,
2023a, this amounts to a reduction by a factor of $\sim$30000\. These terms
are again evaluated using the PySCF module. The prefactors $A$, $\bar{A}$ and
$B$ present in Eq. (16) are straightforwardly implemented in Python according
to Eqs. (17), (18) and (19). Evaluation of the Dyson orbitals needed for the
calculation of the nonadiabatic couplings is implemented as described before
in Ref. Humeniuk  _et al._ , 2013 for arbitrary basis sets for the anion and
the neutral molecule. After construction of the Dyson orbitals from all bound
anionic states to the neutral ground state the nonadiabatic coupling terms are
then calculated according to Eq. (23). To ensure that the wavefunctions of
bound states do not switch their arbitrary signs (which can happen, since the
external quantum-chemistry calculations are independent of each other), the
overlap of electronic wavefunctions of all bound states are tracked throughout
the trajectories and accounted for in all formulae involving the respective
states.

####  III.4.2 Calculation of electronic state coefficients

The electronic degrees of freedom are propagated by solving the time-dependent
Schrödinger equation (5) in the manifold of all considered bound anion and
continuum electronic states using Adams’ method as implemented in the ode
class of Python’s scipy.integrate module Virtanen  _et al._ (2020) with a
user-defined integration time step. For increased computational stability the
equations are beforehand transformed into the interaction picture, introducing
the new electronic state coefficients

$\displaystyle a_{n}(t)=c_{n}(t)\ \mathrm{e}^{\frac{i}{\hbar}H_{nn}t}.$ (41)

Inserting this into Eq. (5) leads to the actually implemented electronic
equation of motion

$\displaystyle\dot{a}_{n}(t)=\sum_{m}\left[-\frac{i}{\hbar}\tilde{H}_{nm}-D_{nm}\right]a_{m}(t)\mathrm{e}^{-\frac{i}{\hbar}(H_{mm}-H_{nn})t}$
(42)

where $\tilde{H}_{nm}$ denotes the Hamiltonian matrix of the system with zeros
on the diagonal.

####  III.4.3 Hopping procedure

Hopping probabilities are directly evaluated according to Eq. (29) from the
state coefficients: A random number between 0 and 1 is generated using the
random function in the numpy.random module and hopping is conducted
accordingly. Once a trajectory hops into a continuum state, it could in
principle be straightforwardly continued on the potential energy surface of
the neutral molecule. Although this can be quite insightful if one is
interested in the subsequent geometric changes of the ionized system, we
follow a different approach and stop the trajectories after electron
detachment since our focus is set on the actual autoionization process. This
allows us to implement a modification of the surface-hopping procedure that
leads to a great improvement of the hopping statistics. The idea is to divide
a single trajectory into ’sub-trajectories’, i.e. to evaluate if a trajectory
hops a number $n_{subtraj}$ of times (see Fig. 1). Every time a sub-trajectory
hops into the continuum, $n_{subtraj}$ is reduced by one and once it reaches
zero, the underlying nuclear dynamics is stopped. It has to be noted that this
procedure is only followed for hops into the continuum, while for hops between
bound anionic states only a single hopping event per trajectory and time step
is possible due to the need to continue the nuclear dynamics on an
unambiguously determined potential energy surface.

###  III.5 Graphical user interface

Figure 2: Left: example page of the graphical input generation tool; right:
output analysis GUI with the electronic state coefficients of a single,
exemplary 2-cyanopyrrolide trajectory. The molecular structure is only an
illustrative image, which was created using the VMD program, and its creation
is not part of the presented program.

Our program package comes with a graphical user interface (GUI) for the input
generation as well as an analysis tool for trajectories. An example of the
former is displayed in Fig. 2. In the input generator, which is started with

⬇

$ hortensia –gui

in addition to all relevant settings for the actual simulation, the user may
find options for the generation of a complete folder structure for the
trajectories as well as bash submit scripts to be used with the Slurm Workload
ManagerSchedMD . Furthermore, the above mentioned Wigner ensemble scripts can
be used and initial conditions can be generated. Therefore it is highly
recommended to use the GUI feature.

Additionally, through the command

⬇

$ hortensia –analysis

one can open the analysis tool which is able to read output files and
visualize them in a sub-window using the matplotlib program package Hunter
(2007).

###  III.6 Installation

The most convenient way to install the program package is downloading or
cloning the repository on our Github pageIssler, Petersen, and Mitrić . In the
main folder, execute

⬇

$ python cysetup.py build_ext –inplace

$ pip install .

to first compile the Cython modules and then install the program. The program
package requires (and will automatically pip install)

  * •

python >= 3.8

  * •

cython \- for faster summation of large arrays, mainly in the calculation of
the two-center integrals in Eqs. (16) and (40)

  * •

scipy \- mainly in the integration of the electronic Schrödinger equation as
outlined in subsection III.4.2

  * •

pyscf \- for the calculation of the two-electron integrals in Eqs. (16) and
(40)

  * •

joblib \- for the parallelization of diabatic couplings

  * •

matplotlib \- for the plots in the sub-window of the analysis tool described
before

and all dependencies thereof. Using the command

⬇

$ pip uninstall hortensia_latest

will uninstall the program package.

##  IV Discussion

In this section we will quantify aspects of the program related to overall
performance. This includes the quality of approximations within the
methodology as well as optimization of time consumption and computational
resources. Moreover the exemplary autoionization dynamics of the
2-cyanopyrrolide anion is discussed.

###  IV.1 Accuracy of k-vector discretization and integral approximations

Figure 3: Comparison of the actual surface area of a unit sphere
($A_{\mathrm{sphere}}=4\pi$, blue line) and the approximated surface area as
described in subsection III.4.1 for up to 104 vector orientations (orange).
The relative error is given in green.

The accuracy of the Fibonacci sphere algorithm for angular discretization in
k-space is illustrated in Fig. 3 by the covered surface area of a unit sphere
using a given number of distributed points. The total surface area (orange
graph) is presented with the relative error
$|A_{\mathrm{fib}}-A_{\mathrm{sphere}}|/A_{\mathrm{sphere}}$ (green graph) to
the exact surface area $4\pi\approx 12.566$ (blue line). The approximated area
rapidly converges to a value of $\sim$12.243, which corresponds to a relative
error of $\sim$2.575 %. Since in the coverage of k-vector lengths no
additional approximation is introduced and for their respective volume
elements the k-space is divided energetically evenly (thus covered exactly
with respect to vector length), the error in the surface area for specific
vector lengths equates to the overall error of the volume elements. Therefore
the sum of these volume elements results in a total volume that deviates by
less than 3 % from the actual sphere for arbitrary numbers of vector
orientations $n_{s}\geq 30$ and lengths $n_{E}$ (giving a total number of wave
vectors $n_{k}=n_{E}\cdot n_{s}$).

Figure 4:  Errors (in %) of hybrid Gaussian-plane wave electron repulsion
integrals $\langle\mathbf{k}_{i}\lambda|\mu\nu\rangle$ for 2-cyanopyrrolide
employing the 6-311++G** +3s2p basis set. The molecular structure has been
optimized in the dipole-bound excited state at the TDDFT/$\omega$B97XD level
using the same basis set. Two types of error measures are reported: a)
$\epsilon_{1}^{\textrm{ap}}=\braket{}{I_{ex}-I_{ap}}{}/\braket{}{I_{ex}}{}$
and b)
$\epsilon_{2}^{\textrm{ap}}~{}=~{}\braket{}{I_{ex}-I_{ap}}{\,/\,|I_{ex}|}^{a}$,
where $I_{ex}$ denotes the exact integral, $I_{ap}$ the approximate value
either according to Eq. (40) (linear, red bars) or assuming the plane wave to
be constant (orange bars). For $\epsilon_{2}^{\textrm{ap}}$ the average has
been computed for all integrals with $|I_{ex}|>10^{-16}\,E_{H}$. To compute
the average the integrals are grouped according to the exponent $\alpha$ of
basis function $\mu$ as core ($\alpha>10\,a_{0}^{-2}$), valence
($10\,a_{0}^{-2}<\alpha<0.1\,a_{0}^{-2}$), and diffuse
($\alpha<0.1\,a_{0}^{-2}$). For each plane wave energy ($E_{1}=0.0015$ eV,
$E_{2}=0.1$ eV, $E_{3}=0.2$ eV), the average has been taken over all distinct
integrals provided by the basis set as well as over 24 different k-vectors
corresponding to the direction vectors of the vertices of a snub cube. In a),
the core and valence error bars are multiplied by a factor of 100 and 10,
respectively, to enhance visibility.

The approximation of the plane wave by the first terms of its Taylor expansion
as introduced in Eq. (40) relies on the assumption that the amplitude of the
plane wave only changes marginally within the extent of the AOs. Fig. 4 shows
a comparison between the approximation with linear terms, an even simpler
constant-wave approximation where
$\textrm{e}^{i\textbf{k}\textbf{r}}\approx\textrm{e}^{i\textbf{k}\textbf{R}_{\mu}}$
and the exact integrals for selected plane wave vectors for 2-cyanopyrrolide,
which serves as an example molecule illustrating the applicability of the
program (see section IV.3 below). Two error measures are compared: a relative
value of average deviations ($\epsilon_{1}$) in Fig. 4a) and an average value
of relative deviations ($\epsilon_{2}$) in Fig. 4b), which differ insofar as
in $\epsilon_{1}$, the deviations between exact and approximate integrals are
averaged first and then divided by the overall average value of the exact
integrals, while in $\epsilon_{2}$, first for each individual integral the
relative error is computed, followed by averaging the results. The averages
are reported for three illustrative plane wave energies and grouped according
to the Gaussian exponent of the basis function sharing its electron coordinate
with the plane wave as "core", "valence" and "diffuse" with decreasing size of
the exponent (for details see Fig. 4). Overall, it becomes evident that for
both error measures, the linear approximation of the plane wave is clearly
superior to the constant approximation. The values of $\epsilon_{1}$ are
always much smaller than those of $\epsilon_{2}$, which is due to the fact
that the relative errors of smaller integrals tend to be larger than those of
larger ones, and the definition of $\epsilon_{1}$ partially compensates for
this fact. Errors larger than a few percent only occur for $\epsilon_{2}$
calculated for diffuse basis functions at larger plane wave energies. Since in
the actual computations, the approximate integrals are employed to calculate
the diabatic couplings and for this, the sum over the entire basis set is
taken (cf. Eq. (16)), especially the smallness of error $\epsilon_{1}$
encourages the use of the linear approximation.

###  IV.2 Optimization of program performance

Where computationally advantageous, we separate the time-dependent and
-independent parts of the underlying equations and pre-calculate the time-
independent terms at the beginning of the simulation. This results in higher
overall memory usage, however of only several hundred MB to a few GB
(depending on the molecular system), but leads to significant time-saving,
which is still a desirable trade-off when calculating on CPU clusters but may
limit the use of the program on single desktop computers.

Furthermore, for increased performance the summation over the four-center
integrals in terms 2 and 4 on the right side of Eq. (16) is implemented as
follows: one first pre-calculates the terms $A$, $\bar{A}$ and $B$ given in
Eqs. (17)-(19) for all AOs. Then the calculation of the four-center integrals
using the PySCF program package is divided into $n_{proc}$ smaller terms,
$n_{proc}$ being the user-defined number of processors, and then evaluated in
parallel utilizing the joblibjob (2021) library by explicit summation over all
AO combinations implemented in a CythonBehnel  _et al._ (2011) module,
therefore reducing the memory usage by ridding oneself from massive arrays
while also improving the runtime performance of this time bottleneck through
parallelization.

Together with the calculation of coupling terms, the most time-consuming step
of the simulation is the two external quantum-chemical calculations needed in
each time step. There are a few options to improve the performance of these
calculations, the easiest of which are to increase the number of utilized
processors and to reduce convergence time by loading the results of the last
time step as an initial guess for the new calculation. Another possibility is
in the choice of basis sets. Finding a basis set for anions prone to
autoionization can be challenging due to the small ionization energies and the
diffusity of the states that comes with it. Therefore one has to consider
basis sets augmented with enough diffuse basis functions to reasonably
describe the properties of the system Herbert (2015). Although popular basis
sets such as doubly and triply augmented Dunning-style basis sets (daug-cc-
pVDZ, taug-cc-pVDZ) are (generally speaking) a potentially good choice for the
description of loosely-bound anions, the size of these basis sets is
computationally prohibitive if one aims to run dynamics simulations and
therefore thousands of consecutive quantum-chemical calculations. A good
alternative can be the usage of smaller basis sets (such as 6-311++G**Clark
_et al._ (1983); Krishnan  _et al._ (1980)) augmented with additional diffuse
functions generated by geometric progressions of the Gaussian exponents as
outlined in Ref. Skurski, Gutowski, and Simons, 2000.

Considering the overall time consumption, no real performance benchmarks exist
with which to compare our program package, since the theory behind it is
rather novel. Therefore we will briefly discuss the specific case of
vinylidene from our work presented in Ref. Issler, Mitrić, and Petersen, 2023a
and the 2-cyanopyrrolide example discussed in detail in section IV.3 below.

The vinylidene dynamics was performed for a total time of 3 ps in 15000
nuclear dynamics time steps at the $\omega$B97XDChai and Head-Gordon
(2008)/d-aug-cc-pVDZ level of theory, which consists of 146 primitive Gaussian
basis functions and 96000 plane waves, amounting to $\sim$460 million
2-electron integrals per time step to be solved (cf. Eqs. (16) and (40)).
Using 6 Intel Xeon E5-2660 (v3) processors per trajectory, the average
computation time was around 11 days and 14 hours with a peak memory usage of
$\sim$9 GB. Of the total time, around 5 days (or 43 %) were needed for the
external quantum-chemistry calculations with the Gaussian09 program package.
It also has to be noted that of the remaining time $\sim$30 % can be
attributed to the calculation of the 2-electron integrals in the diabatic
couplings running on a single processor, which has since been parallelized for
improved performance.

In the simulation of the 2-cyanopyrrolide dynamics, the 6-311++G**+3s2p basis
set consists of 297 primitive Gaussian basis functions which results in
$\sim$7.8 billion 2-electron integrals to be summed over per time step. The
average computation time amounted to 13 days and 12 hours on 10 Intel Xeon
E5-2660 (v3) processors per trajectory, for a total time of 200 fs in 1000
nuclear time steps. The inclusion of an excited state leads to a massive
increase in time consumption in the quantum-chemical calculations (which
account for $\sim$47 %/ 6.4 days of the total computation time) as well as the
evaluation of diabatic couplings, where the summation of all integral terms
(cf. Eq. (16)) is now also conducted on 10 processors using the joblib module.

Table 1:  Comparison of adiabatic electron affinities (AEA), vertical detachment energies to the neutral ground state (VDE${}_{\textrm{GS}}$), vertical attachment energies to the anionic ground state (VAE${}_{\textrm{GS}}$) and excitation energy to the dipole-bound state ($\Delta$E${}_{\textrm{DBS}}$) for the 2-cyanopyrrolide anion. The superscript indicates at which optimized geometry (a = anion, n = neutral) the respective value is calculated. The method used in the dynamics simulation is indicated in bold font. The added basis functions +Xs etc. are generated according to Ref. Skurski, Gutowski, and Simons, 2000 (with a factor for the geometric progression of 3.5) and centered on the nitrogen atoms. All energies values are given in eV. In some cases AEA${}_{\textrm{DBS}}$ could not be obtained, instead the vertical attachment energy at the neutral equilibrium geometry is given as an approximation (denoted with *). ${}^{\textrm{a}}$ Ref. Dunning, 1989, ${}^{\textrm{b}}$ Ref. Dunning, 1989; Kendall, Dunning, and Harrison, 1992, ${}^{\textrm{c}}$ Ref. Dunning, 1989; Kendall, Dunning, and Harrison, 1992; Woon and Dunning, 1994, ${}^{\textrm{d}}$ Ref. Clark  _et al._ , 1983, ${}^{\textrm{e}}$ Ref. Clark  _et al._ , 1983; Krishnan  _et al._ , 1980, ${}^{\textrm{f}}$ Ref. Purvis and Bartlett, 1982; Scuseria, Janssen, and Schaefer, 1988 Method | AEA${}_{\textrm{GS}}$ | AEA${}_{\textrm{DBS}}$ | VDE${}^{\textrm{a}}_{\textrm{GS}}$ | VAE${}^{\textrm{n}}_{\textrm{GS}}$ |  $\Delta$E${}^{\textrm{n}}_{\textrm{DBS}}$
---|---|---|---|---|---
$\omega$B97XD / aug-cc-pVDZ${}^{\textrm{b}}$ | 3.075 | -0.674 | 3.225 | 2.932 | 3.612
$\omega$B97XD / d-aug-cc-pVDZ${}^{\textrm{c}}$ | 3.071 | -0.117* | 3.221 | 2.929 | 3.046
$\omega$B97XD / t-aug-cc-pVDZ${}^{\textrm{c}}$ | 3.070 | 0.031* | 3.220 | 2.928 | 2.897
$\omega$B97XD / t-aug-cc-pVTZ${}^{\textrm{c}}$ | 3.044 | 0.044* | 3.200 | 2.899 | 2.855
$\omega$B97XD / 6-31+G**${}^{\textrm{d}}$ | 3.062 | -1.374 | 3.212 | 2.919 | 4.302
$\omega$B97XD / 6-31+G** +3s2p | 3.064 | 0.060* | 3.217 | 2.921 | 2.861
$\omega$B97XD / 6-311++G**${}^{\textrm{e}}$ | 3.095 | -0.759 | 3.249 | 2.949 | 3.716
$\omega$B97XD / 6-311++G** +3s | 3.095 | 0.064 | 3.249 | 2.949 | 2.887
$\omega$B97XD / 6-311++G** +3s2p | 3.094 | 0.063 | 3.248 | 2.948 | 2.887
$\omega$B97XD / 6-311++G** +3s2p2d | 3.076 | 0.055* | 3.230 | 2.930 | 2.875
EOM-CCSD${}^{\textrm{f}}$ / t-aug-cc-pVDZ | 2.844 | -0.102* | 3.016 | 2.686 | 2.788
EOM-CCSD / 6-31+G** +3s2p | 2.679 | -0.346* | 2.855 | 2.518 | 2.864
EOM-CCSD / 6-311++G** +3s2p | 2.748 | -0.148* | 2.929 | 2.584 | 2.732
ExperimentYuan, Zhang, and Wang (2022) | 3.0981 | 0.0298 |  |  |
Figure 5: HONTO and LUNTO of 2-cyanopyrrolide at the optimized geometry of the
dipole-bound first excited state at the $\omega$B97XD/ 6-311++G** + 3s2p level
of theory with an isovalue of 0.003. Figure 6: a) Overlay of all initial
molecular structures used in the dynamics simulation of 2-cyanopyrrolide, b)
distribution for 10000 initial conditions as a function of distance (in Å)
between the nitrogen (1) and carbon (3) atom as marked in a), showing a
bimodal structure. Figure 7: a) Mollweide projection of the angular
distribution of ejected electrons in the 2-cyanopyrrolide dynamics, summed
over all energies. The x-axis ($\varphi=0$, $\theta=90$ degrees) is aligned
with the cyano group and the molecule lies within the xy-plane ($\theta=90$
degrees); b) Slices through the Mollweide projection at $\varphi$ angles of 0
(positive x direction, blue), 180 (negative x direction, orange), 90 (positive
y direction, green) and 270 (negative y direction, red) degrees. Figure 8: a)
Simulated electron kinetic energy distribution of all hopping events after
excitation of mode $\nu_{11}$ and propagation for 200 fs (orange histogram)
and running average over 5 points/7.5 meV (red curve), b) time-dependent
population of all bound anion states (dark green) and exponential fit with a
time constant of $\tau=500$ fs (light green).

###  IV.3 Illustrative example: Autoionization of the 2-cyanopyrrolide anion

To illustrate the scope of our program, we simulated the vibration-induced
autoionization dynamics of the example anion 2-cyanopyrrolide. Experimentally,
this molecule was measured to have an adiabatic electron affinity of 3.0981 eV
and possesses a Rydberg-s type dipole-bound state 29.8 meV below the
ionization threshold.Yuan, Zhang, and Wang (2022) As can be seen in Table 1,
which compares several quantum-chemistry methods and basis sets with the data
measured by Wang $et\ al.$, the experimental data is reproduced quite well
using the $\omega$B97XD functional and large, diffuse basis sets such as
triply augmented pVDZ/pVTZ. Moreover, although the description of the molecule
with standard Pople-type basis sets is fairly inaccurate, further augmentation
with extra diffuse basis functions (see Ref. Skurski, Gutowski, and Simons,
2000), in this case placed on the nitrogen atoms, also leads to good agreement
with the experimental values. At the same time this approach retains a
significantly smaller total number of basis functions, therefore keeping
computational effort manageable. Fig. 5 shows the HONTO and LUNTO, visualizing
the spatial distribution of the excess electron in the ground and excited
state at the optimized geometry of the dipole-bound first excited state
employing the $\omega$B97XD functional and the 6-311++G** basis set augmented
with three diffuse s- and two diffuse p-functions on each nitrogen atom
(henceforth abbreviated as 6-311++G** + 3s2p). The shape of the excess
electron’s probability distribution in the dipole-bound state is of s-type,
showing that employing additional higher polarization functions (d-/f-type)
would lead to no further improvement in the description of the system. This is
in complete agreement with a dipole moment of the neutral species of 5.02 D,
well below the second critical dipole moment of $\sim$10 D needed for the
binding of an electron in a p-type orbital,Jordan and Wang (2003) consequently
resulting in an s-type distribution centered around the positive end of the
molecular dipole vector.

Using the 6-311++G** + 3s2p basis set with the $\omega$B97XD functional, we
simulated the vibration-induced autoionization dynamics in the first excited
state with the normal mode at 946 cm-1 of A’ symmetry ($\nu_{11}$ when sorted
by increasing mode energy irrespective of symmetry) excited by one vibrational
quantum. The initial conditions were generated as described in subsection
III.2. Mode $\nu_{11}$ involves a symmetric stretching of the C-H bonds at
carbon atoms 4 and 5 as well as a ring breathing motion affecting mostly the
ring N and carbon 3. The numbering of atoms is provided in Fig. 6a), which
illustrates the resulting set of initial conditions by the superposition of
all initial structures. In Fig. 6b), the distance between the ring nitrogen
and the carbon 3 is depicted, which exhibits a bimodal distribution typical
for an excited vibrational state. The particular choice of vibrational
excitation corresponds to the experimentally observed resonance 7 of the
photodetachment spectrum in Ref. Yuan, Zhang, and Wang, 2022. The simulation
was carried out propagating an ensemble of 53 trajectories for a total of 200
fs (1000 nuclear time steps) with a discretized continuum of 400 plane wave
energies evenly spaced from 0.0 eV to 0.138 eV and 96 orientations per energy.
The maximum allowed kinetic energy of the plane wave is the sum of the
vibrational excitation energy and the difference in zero-point energies of
anion and neutral system, that is, the maximum excess energy available upon
ionization.

Notice that due to the very low electron binding energy of the dipole-bound
state and the approximative nature of the quantum chemically determined
energies, it is challenging to precisely reproduce subtle binding energy
differences on the meV scale along the trajectories. Thus, some instances of
negative VDE occur in the dynamics. However, the experimental data from Ref.
Yuan, Zhang, and Wang, 2022 only feature a peak attributed to vibrational
autoionization. Therefore, we only include the latter in our simulation and
neglect adiabatic ionization.

The nuclear dynamics following the vibrational excitation is characterized by
relatively small amplitude motion. This is due to the overall low internal
energy of the molecule and its rigidity as a cyclic system. In the course of
the dynamics, the molecular dipole moment associated with the neutral core,
which is responsible for electron binding in the excited state, exhibits
slight oscillatory behavior while being approximately situated in the
molecular plane. This leads to an anisotropic ejection of electrons
predominantly in the molecular plane along the axis containing the cyano
group, as can be inferred from the Mollweide projection of the angle-dependent
distribution of k-vectors, summed over all k values shown in Fig. 7a). The
resulting electron distribution is thus p-shaped, with maxima along the
x-(cyano group) axis and minima in the yz-plane exhibiting only about 20% of
the maximal intensity, as can be seen in Fig. 7b). This observation is in line
with the qualitative considerations of nonadiabatic autoionization from
dipole-bound states outlined in Ref. Simons, 2020. No transitions to the
anionic ground state are observed in our simulation due to a large energy gap
regardless of geometry, therefore the angular electron distribution is solely
due to ionization from the s-type dipole-bound state.

Regarding the electron kinetic energies, the distribution displayed in Fig.
8a) is obtained, exhibiting a broad peak near the maximally possible energy of
0.138 eV. This can be attributed to a transition in which the vibrational
energy of the excited mode is transferred completely to the outgoing electron,
i.e., the vibrational energy of the molecule is reduced by one quantum in line
with the propensity rules for vibrational autoionization established by
SimonsSimons (1981). Further analysis of the peak shape should be taken with
care, since for conceptual reasons vibrational resolution is not within the
scope of quantum-classical dynamics.

Besides the spatial and energetic distribution of the ejected electrons, our
simulation provides access to the timescale in which the ionization process
takes place. Fig. 8b) shows the time-dependent population of the bound anionic
states, which exhibits a rapid decay that can be fit to an exponential
function with a time constant of 500 fs. This value corresponds to a spectral
width of around 70 cm-1, which is of comparable size to the observation made
in Ref. Yuan, Zhang, and Wang, 2022.

Overall this example calculation shows the applicability and scope of the
method in the context of small to medium sized molecular anions, providing a
means to gain molecular-level insight into the spatio-temporal dynamics of
vibration-induced autoionization processes complementary to experimental
measurements.

##  V Conclusion

We have presented the Python program package HORTENSIA (Hopping real-time
trajectories for electron-ejection by nonadiabatic self-ionization in anions)
for the simulation of vibration-induced autoionization processes in molecular
anions. The program implements our recently introduced extended surface
hopping approach for the quantum-classical description of nonadiabatic
autoionization dynamics, where the electronic degrees of freedom are treated
quantum-mechanically, while the nuclear motion is represented by classical
trajectories. The electronic states included in the dynamics simulation
comprise the bound adiabatic anionic states and discretized ’ionized system’
states composed of a neutral core and a free electron wave function, between
which nonadiabatic transitions are simulated in a stochastical manner from
hopping probabilities obtained from changes in electronic state coefficients
according to Tully’s fewest-switches algorithm. The time-dependent state
coefficients are calculated by solution of the electronic Schrödinger equation
containing the nonadiabatic as well as diabatic couplings between the
considered electronic states according to our presented methodology.

As shown in the example of 2-cyanopyrrolide, time- and angle-resolved electron
kinetic energy signals are obtained directly from the surface-hopping
trajectories. Since no deactivation to the ground state is observed in our
simulation, autoionization with a time constant of 500 fs is identified as the
only available deactivation pathway in the dipole-bound state of
2-cyanopyrrolide on the simulated timescale, with an anisotropic, p-like
ejection of electrons along the cyano-axis. Moreover, with our program
geometric data is yielded which allows for the structural analysis of
molecules throughout the autoionization dynamics, providing easy access to
geometric characteristics of the considered system, as demonstrated
extensively in the example of the vinylideneIssler, Mitrić, and Petersen
(2023a) and 1-nitropropaneIssler, Mitrić, and Petersen (2023b) anions.

Furthermore, the implementation and internal structure of our program package
was discussed, which also consists of secondary functionalities such as an
input generator and a routine for the creation of initial conditions for
nuclear coordinates and velocities within an easy-to-operate graphical user
interface (GUI). Moreover, the program package provides the user with an
additional GUI for the analysis and graphical representation of the most
important dynamics results.

In the future, useful extensions of the methodology could be the
implementation of neutral molecules to be ionized, which requires the
description of scattering states interacting with a cationic core, as well as
the inclusion of laser field coupling (analogous to the FISH methodMitrić,
Petersen, and Bonačić-Koutecký (2009)) to describe photoionization beyond the
perturbative limit, thereby providing an extension of the approach developed
in Ref. Humeniuk  _et al._ , 2013. In addition, the treatment of
electronically adiabatic autoionization could be combined with an ab inito
computation of the electronic resonance lifetimes, e.g., along the lines
presented in Ref. Gyamfi and Jagau, 2022.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Author declarations

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Kevin Issler: Data curation (lead); Formal analysis (lead); Investigation
(lead); Methodology (equal); Software (lead); Visualization (lead); Writing -
original draft (equal). Roland Mitrić: Conceptualization (equal); Funding
acquisition (lead); Methodology (supporting); Project administration (lead);
Resources (lead); Supervision (equal); Writing - review & editing (equal).
Jens Petersen: Conceptualization (equal); Formal analysis (supporting);
Methodology (lead); Software (supporting); Supervision (equal); Visualization
(supporting); Writing - original draft (equal), Writing - review & editing
(equal).

## References

  * Illenberger and Momigny (1992) E. Illenberger and J. Momigny,  _Gaseous Molecular Ions. An Introduction to Elementary Processes Induced by Ionization_ (Steinkopff Verlag, Darmstadt; Springer-Verlag, New York, 1992).
  * Boudaïffa  _et al._ (2000) B. Boudaïffa, P. Cloutier, D. Hunting, M. A. Huels, and L. Sanche, “Resonant formation of DNA strand breaks by low-energy (3 to 20 eV) electrons,” Science 287, 1658 (2000).
  * Martin  _et al._ (2004) F. Martin, P. D. Burrow, Z. Cai, P. Cloutier, D. Hunting, and L. Sanche, “DNA strand breaks induced by 0 - 4 eV electrons: The role of shape resonances,” Phys. Rev. Lett. 93, 068101 (2004).
  * Arumainayagam  _et al._ (2010) C. R. Arumainayagam, H.-L. Lee, R. B. Nelson, D. R. Haines, and R. P. Gunawardane, “Low-energy electron-induced reactions in condensed matter,” Surf. Sci. Rep. 65, 1–44 (2010).
  * Liu  _et al._ (2013) H.-T. Liu, C.-G. Ning, D.-L. Huang, P. D. Dau, and L.-S. Wang, “Observation of Mode-Specific Vibrational Autodetachment from Dipole-Bound States of Cold Anions,” Angew. Chem. Int. Ed. 52, 8976–8979 (2013).
  * Liu  _et al._ (2014) H.-T. Liu, C.-G. Ning, D.-L. Huang, and L.-S. Wang, “Vibrational Spectroscopy of the Dehydrogenated Uracil Radical by Autodetachment of Dipole-Bound Excited States of Cold Anions,” Angew. Chem. Int. Ed. 53, 2464–2468 (2014).
  * Zhu, Liu, and Wang (2017) G.-Z. Zhu, Y. Liu, and L.-S. Wang, “Observation of Excited Quadrupole-Bound States in Cold Anions,” Phys. Rev. Lett. 119, 023002 (2017).
  * Yuan, Zhang, and Wang (2022) D.-F. Yuan, Y.-R. Zhang, and L.-S. Wang, “Dipole-Bound State, Photodetachment Spectroscopy, and Resonant Photoelectron Imaging of Cryogenically-Cooled 2-Cyanopyrrolide,” J. Phys. Chem. A 126, 6416–6428 (2022).
  * Harbach  _et al._ (2013) P. H. P. Harbach, M. Schneider, S. Faraji, and A. Dreuw, “Intermolecular coulombic decay in biology: The initial electron detachment from FADH- in DNA photolyases,” J. Phys. Chem. Lett. 4, 943–949 (2013).
  * Adams, Schneider, and Weber (2010) C. L. Adams, H. Schneider, and J. M. Weber, “Vibrational autodetachment - Intramolecular vibrational relaxation translated into electronic motion,” J. Phys. Chem. A 114, 4017–4030 (2010).
  * Bull, West, and Verlet (2016) J. N. Bull, C. W. West, and J. R. R. Verlet, “Ultrafast dynamics of formation and autodetachment of a dipole-bound state in an open-shell $\pi$-stacked dimer anion,” Chem. Sci. 7, 5352–5361 (2016).
  * DeVine  _et al._ (2017) J. A. DeVine, M. L. Weichman, B. Laws, J. Chang, M. C. Babin, G. Balerdi, C. Xie, C. L. Malbon, W. C. Lineberger, D. R. Yarkony, R. W. Field, S. T. Gibson, J. Ma, H. Guo, and D. M. Neumark, “Encoding of vinylidene isomerization in its anion photoelectron spectrum,” Science 358, 336 (2017).
  * DeVine  _et al._ (2018) J. A. DeVine, M. L. Weichman, C. Xie, M. C. Babin, M. A. Johnson, J. Ma, H. Guo, and D. M. Neumark, “Autodetachment from vibrationally excited vinylidene anions,” J. Phys. Chem. Lett. 9, 1058 (2018).
  * Adams, Hansen, and Weber (2019) C. L. Adams, K. Hansen, and J. M. Weber, “Vibrational Autodetachment from Anionic Nitroalkane Chains: From Molecular Signatures to Thermionic Emission,” J. Phys. Chem. A 123, 8562–8570 (2019).
  * Anstöter  _et al._ (2020) C. S. Anstöter, G. Mensa-Bonsu, P. Nag, M. Ranković, R. Kumar T. P., A. N. Boichenko, A. V. Bochenkova, J. Fedor, and J. R. R. Verlet, “Mode-specific vibrational autodetachment following excitation of electronic resonances by electrons and photons,” Phys. Rev. Lett. 124, 203401 (2020).
  * Naff, Cooper, and Compton (1968) W. T. Naff, C. D. Cooper, and R. N. Compton, “Transient negative-ion states in alicyclic and aromatic fluorocarbon molecules,” J. Chem. Phys. 49, 2784–2788 (1968).
  * Jordan and Burrow (1987) K. D. Jordan and P. D. Burrow, “Temporary anion states of polyatomic hydrocarbons,” Chem. Rev. 87, 557–588 (1987).
  * Suess, Parthasarathy, and Dunning (2002) L. Suess, R. Parthasarathy, and F. B. Dunning, “Nondissociative low-energy electron attachment to SF6, C6F6, C10F8, and c-C7F14: Negative ion lifetimes,” J. Chem. Phys. 117, 11222–11227 (2002).
  * Bull, Anstöter, and Verlet (2019) J. N. Bull, C. S. Anstöter, and J. R. R. Verlet, “Ultrafast valence to non-valence excited state dynamics in a common anionic chromophore,” Nat. Commun. 10, 5820 (2019).
  * Verlet  _et al._ (2020) J. R. R. Verlet, C. S. Anstöter, J. N. Bull, and J. P. Rogers, “Role of Nonvalence States in the Ultrafast Dynamics of Isolated Anions,” J. Phys. Chem. A 124, 3507–3519 (2020).
  * Kang, An, and Kim (2020) D. H. Kang, S. An, and S. K. Kim, “Real-Time Autodetachment Dynamics of Vibrational Feshbach Resonances in a Dipole-Bound State,” Phys. Rev. Lett. 125, 093001 (2020).
  * Kang  _et al._ (2021) D. H. Kang, J. Kim, M. Cheng, and S. K. Kim, “Mode-specific autodetachment dynamics of an excited non-valence quadrupole-bound state,” J. Phys. Chem. Lett. 12, 1947–1954 (2021).
  * Herbert (2015) J. M. Herbert, “The quantum chemistry of loosely-bound electrons. in reviews in computational chemistry volume 28,” (John Wiley & Sons, Ltd, 2015) Chap. 8, pp. 391–517.
  * Berry (1966) R. S. Berry, “Ionization of molecules at low energies,” J. Chem. Phys. 45, 1228–1245 (1966).
  * Acharya, Kendall, and Simons (1984) P. K. Acharya, R. A. Kendall, and J. Simons, “Vibration-induced electron detachment in molecular anions,” J. Am. Chem. Soc. 106, 3402–3407 (1984).
  * Simons (1999) J. Simons, “Time-domain and tunneling pictures of nonadiabatic induced electron ejection in molecular anions,” J. Phys. Chem. A 103, 9408–9416 (1999).
  * Simons (1981) J. Simons, “Propensity rules for vibration-induced electron detachment of anions,” J. Am. Chem. Soc. 103, 3971 (1981).
  * Humeniuk  _et al._ (2013) A. Humeniuk, W. Wohlgemuth, T. Suzuki, and R. Mitrić, “Time-resolved photoelectron imaging spectra from non-adiabatic molecular dynamics simulations,” J. Chem. Phys. 139, 134104 (2013).
  * Persico and Granucci (2014) M. Persico and G. Granucci, “An overview of nonadiabatic dynamics simulations methods, with focus on the direct approach versus the fitting of potential energy surfaces,” Theor. Chem. Acc. 133, 1526 (2014).
  * Richings  _et al._ (2015) G. Richings, I. Polyak, K. Spinlove, G. Worth, I. Burghardt, and B. Lasorne, “Quantum dynamics simulations using gaussian wavepackets: the vMCG method,” Int. Rev. Phys. Chem. 34, 269–308 (2015).
  * Crespo-Otero and Barbatti (2018) R. Crespo-Otero and M. Barbatti, “Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics,” Chem. Rev. 118, 7026–7068 (2018).
  * Röder  _et al._ (2019) A. Röder, J. Petersen, K. Issler, I. Fischer, R. Mitrić, and L. Poisson, “Exploring the excited-state dynamics of hydrocarbon radicals, biradicals and carbenes using time-resolved photoelectron spectroscopy and field-induced surface hopping simulations,” J. Phys. Chem. A 123, 10643 (2019).
  * Barbatti  _et al._ (2014) M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico, and H. Lischka, “Newton-X: a surface-hopping program for nonadiabatic molecular dynamics,” WIREs Comput. Mol. Sci. 4, 26–33 (2014).
  * Mai, Marquetand, and González (2018) S. Mai, P. Marquetand, and L. González, “Nonadiabatic dynamics: The sharc approach,” WIREs Comp. Mol. Sci. 8, e1370 (2018).
  * Du and Lan (2015) L. Du and Z. Lan, “An on-the-fly surface-hopping program jade for nonadiabatic molecular dynamics of polyatomic systems: Implementation and applications,” J. Chem. Theory Comput. 11, 1360–1374 (2015).
  * Issler, Mitrić, and Petersen (2023a) K. Issler, R. Mitrić, and J. Petersen, “Quantum-classical Dynamics of Vibration-Induced Autoionization in Molecules,” J. Chem. Phys. 158, 034107 (2023a).
  * Tully (1990) J. C. Tully, “Molecular dynamics with electronic transitions,” J. Chem. Phys. 93, 1061–1071 (1990).
  * Jordan and Wang (2003) K. D. Jordan and F. Wang, “Theory of dipole-bound anions,” Annu. Rev. Phys. Chem. 54, 367–396 (2003).
  * Issler, Mitrić, and Petersen (2023b) K. Issler, R. Mitrić, and J. Petersen, “A trajectory surface hopping study of the vibration-induced autoionization dynamics of the 1-nitropropane anion,” Theo. Chem. Acc. , submitted (2023b).
  * Mitrić, Werner, and Bonačić-Koutecký (2008) R. Mitrić, U. Werner, and V. Bonačić-Koutecký, “Nonadiabatic dynamics and simulation of time resolved photoelectron spectra within time-dependent density functional theory: Ultrafast photoswitching in benzylideneaniline,” J. Chem. Phys. 129, 164118 (2008).
  * Werner  _et al._ (2008) U. Werner, R. Mitrić, T. Suzuki, and V. Bonačić-Koutecký, “Nonadiabatic dynamics within the time dependent density functional theory: Ultrafast photodynamics in pyrazine,” Chem. Phys. 349, 319–324 (2008).
  * Werner, Mitrić, and Bonačić-Koutecký (2010) U. Werner, R. Mitrić, and V. Bonačić-Koutecký, “Simulation of time resolved photoelectron spectra with stieltjes imaging illustrated on ultrafast internal conversion in pyrazine,” J. Chem. Phys. 132, 174301 (2010).
  * Martin (2003) R. L. Martin, “Natural transition orbitals,” J. Chem. Phys. 118, 4775–4777 (2003).
  * McMurchie and Davidson (1978) L. E. McMurchie and E. R. Davidson, “One- and two-electron integrals over cartesian Gaussian functions,” J. Comp. Phys. 26, 218 (1978).
  * Mitrić, Petersen, and Bonačić-Koutecký (2011) R. Mitrić, J. Petersen, and V. Bonačić-Koutecký, “Multistate nonadiabatic dynamics “on the fly” in complex systems and its control by laser fields,” in  _Conical Intersections - Theory, Computation and Experiment_, Advanced Series in Physical Chemistry, Vol. 17 (World Scientific, Singapore, 2011) p. 497–568.
  * Frisch  _et al._ (2016) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, “Gaussian 16 Revision A.03,”  (2016), Gaussian Inc. Wallingford CT.
  * Shao  _et al._ (2015) Y. Shao, Z. Gan, E. Epifanovsky, A. T. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kuś, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. W. III, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. D. Jr., H. Do, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. Hanson-Heine, P. H. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. D. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.-P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, E. Neuscamman, C. M. Oana, R. Olivares-Amaya, D. P. O’Neill, J. A. Parkhill, T. M. Perrine, R. Peverati, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, S. M. Sharada, S. Sharma, D. W. Small, A. Sodt, T. Stein, D. Stück, Y.-C. Su, A. J. Thom, T. Tsuchimochi, V. Vanovschi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, J. Yang, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhao, B. R. Brooks, G. K. Chan, D. M. Chipman, C. J. Cramer, W. A. G. III, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. S. III, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. V. Voorhis, J. M. Herbert, A. I. Krylov, P. M. Gill, and M. Head-Gordon, “Advances in molecular quantum chemistry contained in the Q-Chem 4 program package,” Mol. Phys. 113, 184–215 (2015).
  * Skurski, Gutowski, and Simons (2000) P. Skurski, M. Gutowski, and J. Simons, “How to choose a one-electron basis set to reliably describe a dipole-bound anion,” Int. J. Quantum Chem. 80, 1024–1038 (2000).
  * Swope  _et al._ (1982) W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, “A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters,” J. Chem. Phys. 76, 637 (1982).
  * Humphrey, Dalke, and Schulten (1996) W. Humphrey, A. Dalke, and K. Schulten, “VMD – Visual Molecular Dynamics,” J. Mol. Graph. 14, 33–38 (1996).
  * Swinbank and James Purser (2006) R. Swinbank and R. James Purser, “Fibonacci grids: A novel approach to global modelling,” Q. J. R. Meteorol. Soc. 132, 1769–1793 (2006).
  * Sun (2015) Q. Sun, “Libcint: An efficient general integral library for gaussian basis functions,” J. Comput. Chem. 36, 1664–1671 (2015).
  * Sun  _et al._ (2018) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, “PySCF: the Python-based simulations of chemistry framework,” WIREs Comput. Mol. Sci. 8, e1340 (2018).
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, “Recent developments in the PySCF program package,” J. Chem. Phys. 153, 024109 (2020).
  * Colle, Fortunelli, and Simonucci (1987) R. Colle, A. Fortunelli, and S. Simonucci, “A mixed basis set of plane waves and Hermite Gaussian functions. Analytic expressions of prototype integrals,” Il Nuovo Cimento D 9, 969–977 (1987).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nat. Methods 17, 261–272 (2020).
  * (57) SchedMD, “Slurm Workload Manager,” https://slurm.schedmd.com/.
  * Hunter (2007) J. D. Hunter, “Matplotlib: A 2D graphics environment,” Comput. Sci. Eng. 9, 90–95 (2007).
  * (59) K. Issler, J. Petersen, and R. Mitrić, “HORTENSIA program package,” https://github.com/mitric-lab/HORTENSIA_LATEST.git.
  * job (2021) “Joblib: running Python functions as pipeline jobs. Release 1.1.0,”  (2021).
  * Behnel  _et al._ (2011) S. Behnel, R. Bradshaw, C. Citro, L. Dalcin, D. S. Seljebotn, and K. Smith, “Cython: The best of both worlds,” Comput. Sci. Eng. 13, 31–39 (2011).
  * Clark  _et al._ (1983) T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. V. R. Schleyer, “Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li–F,” J. Comput. Chem. 4, 294 (1983).
  * Krishnan  _et al._ (1980) R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, “Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions,” J. Chem. Phys. 72, 650 (1980).
  * Chai and Head-Gordon (2008) J.-D. Chai and M. Head-Gordon, “Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections,” Phys. Chem. Chem. Phys. 10, 6615 (2008).
  * Dunning (1989) T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen,” J. Chem. Phys. 90, 1007 (1989).
  * Kendall, Dunning, and Harrison (1992) R. A. Kendall, T. H. Dunning, and R. J. Harrison, “Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions,” J. Chem. Phys. 96, 6796 (1992).
  * Woon and Dunning (1994) D. E. Woon and T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties,” J. Chem. Phys. 100, 2975 (1994).
  * Purvis and Bartlett (1982) G. D. Purvis and R. J. Bartlett, “A full coupled-cluster singles and doubles model: The inclusion of disconnected triples,” J. Chem. Phys. 76, 1910 (1982).
  * Scuseria, Janssen, and Schaefer (1988) G. E. Scuseria, C. L. Janssen, and H. F. Schaefer, “An efficient reformulation of the closed-shell coupled cluster single and double excitation (CCSD) equations,” J. Chem. Phys. 89, 7382 (1988).
  * Simons (2020) J. Simons, “Ejecting electrons from molecular anions via shine, shake/rattle, and roll,” J. Phys. Chem. A 124, 8778–8797 (2020).
  * Mitrić, Petersen, and Bonačić-Koutecký (2009) R. Mitrić, J. Petersen, and V. Bonačić-Koutecký, “Laser-field-induced surface-hopping method for the simulation and control of ultrafast photodynamics,” Phys. Rev. A 79, 053416 (2009).
  * Gyamfi and Jagau (2022) J. A. Gyamfi and T.-C. Jagau, “Ab initio molecular dynamics of temporary anions using complex absorbing potentials,” J. Phys. Chem. Lett. 13, 8477–8483 (2022).