# Self-Adaptive Real-Time Time-Dependent Density Functional Theory for X-ray
Absorptions

Linfeng Ye  Qingdao Institute for Theoretical and Computational Sciences,
Institute of Frontier and Interdisciplinary Science, Shandong University,
Qingdao, Shandong 266237, P. R. China  Hao Wang<EMAIL_ADDRESS>Qingdao
Institute for Theoretical and Computational Sciences, Institute of Frontier
and Interdisciplinary Science, Shandong University, Qingdao, Shandong 266237,
P. R. China  Yong Zhang  Qingdao Institute for Theoretical and Computational
Sciences, Institute of Frontier and Interdisciplinary Science, Shandong
University, Qingdao, Shandong 266237, P. R. China  Wenjian Liu
<EMAIL_ADDRESS>Qingdao Institute for Theoretical and Computational Sciences,
Institute of Frontier and Interdisciplinary Science, Shandong University,
Qingdao, Shandong 266237, P. R. China

###### Abstract

Real-time time-dependent density functional theory (RT-TDDFT) can in principle
access the whole absorption spectrum of a many-electron system exposed to a
narrow pulse. However, this requires an accurate and efficient propagator for
the numerical integration of the time-dependent Kohn-Sham equation. While a
low-order time propagator is already sufficient for the low-lying valence
absorption spectra, it is no longer the case for the X-ray absorption spectra
(XAS) of systems composed even only of light elements, for which the use of a
high-order propagator is indispensable. It is then crucial to choose a largest
possible time step and a shortest possible simulation time, so as to minimize
the computational cost. To this end, we propose here a robust AutoPST approach
to determine automatically (Auto) the propagator (P), step (S), and time (T)
for relativistic RT-TDDFT simulations of XAS.

††preprint: AIP/123-QED

##  I Introduction

Triggered by the development of advanced synchrotron light sources and X-ray
free electron lasers, X-ray absorption spectroscopy (XAS) has become an
increasingly important technique for probing local electronic and geometric
structure of matter. In XAS, innermost core electrons are excited to bound
valence or continuum states in energetically distinct absorption edges, which
are conventionally labeled according to the origins of the electronic
transitions ($K$-edge for $1s$, $L_{1}$-edge for $2s$, $L_{2}$-edge for
$2p_{1/2}$, and $L_{3}$-edge for $2p_{3/2}$, etc). There are two main regions
in an XAS spectrum, X-ray absorption near-edge structure (XANES) near the
rising edge with sharp resonance peaks and extended X-ray absorption fine
structure (EXAFS) after the XANES region with gentle oscillations. The former
provides electronic fingerprints of a particular atom (e.g., oxidation state,
coordination, and bonding, etc), whereas the latter gives structural
information about neighboring atoms (e.g., identity, distances, and
coordination/solvation shells, etc). While the essential quantities
(scattering amplitudes and phase shifts) in EXAFS can well be modeled already
by a simple, damped spherical photoelectron wave function approximation, the
complicated features of XANES spectra require more sophisticated
treatmentsNorman and Dreuw (2018), among which both static and time-dependent
density functional theory (TDDFT) are widely adopted due to their good
tradeoff between accuracy and efficiency. The former is commonly referred to
as delta self-consistent field ($\Delta$SCF) and has a very appealing feature:
orbital relaxation (which is of vast importance for core excitations) is fully
accounted for by SCF calculation of each core excited configuration. In
particular, full spin adaptation can be achieved within the framework of
multi-state DFTGao  _et al._ (2016); Grofe  _et al._ (2017); Yang  _et al._
(2019); Zhao  _et al._ (2021). Yet, such one-by-one calculations of core
excited states may not always be possible due to the non-Aufbau nature, given
the availability of many algorithmsGilbert  _et al._ (2008); Liu  _et al._
(2014); Derricotte and Evangelista (2015); Liang  _et al._ (2017); Derricotte
and Evangelista (2017); Oakley and Klobukowski (2018); Ehlert and Klamroth
(2020); Carter-Fenk and Herbert (2020); Shea  _et al._ (2020); Hait and Head-
Gordon (2020); Grofe  _et al._ (2020). In contrast, TDDFT can access all core
excited states in one shot, by either diagonalization as in linear response
(LR) TDDFTFronzoni  _et al._ (2004); Stener  _et al._ (2003); Ray  _et al._
(2007); Besley and Asmuruf (2010); Liang  _et al._ (2011); Van Kuiken  _et
al._ (2012); Li and Liu (2010); Li  _et al._ (2011); Li and Liu (2011); Van
Kuiken  _et al._ (2013); Verma  _et al._ (2016); Peng  _et al._ (2019); Nenov
_et al._ (2019) or spectral analysis of the time signal generated by real-time
(RT) TDDFT Lopata  _et al._ (2012); Lopata and Govind (2013); Fernando  _et
al._ (2015); Sissay  _et al._ (2016); Bruner  _et al._ (2016); Neville and
Schuurman (2018); Chen and Lopata (2020); Li  _et al._ (2020); Yang  _et al._
(2022). Yet, both LR-TDDFT and RT-TDDFT lacks orbital relaxation, thereby
yielding XANES spectra that usually have large shifts relative to the
experimental ones. Since core excitations look very much like charge transfer
excitations, it is clear that such large shifts can be alleviated to a large
extent by using tuned range-separated functionalsFernando  _et al._ (2015);
Verma and Bartlett (2016); Zapata  _et al._ (2019). Anyhow, peak separations
predicted by LR-TDDFT and RT-TDDFT are usually very good. As such, such one-
for-all approaches are much more appealing than the one-by-one $\Delta$SCF
approaches, at least from the computational point of view. Naturally, the
relativistic counterparts of LR-TDDFTGao  _et al._ (2004, 2005); Bast  _et
al._ (2009); Komorovsky  _et al._ (2019); Konecny  _et al._ (2019); Wang  _et
al._ (2005); Peng  _et al._ (2005); Xu  _et al._ (2009a, b); Kühn and Weigend
(2013, 2014); Mori  _et al._ (2014); Egidi  _et al._ (2016, 2017); Stetina
_et al._ (2019); Liu and Xiao (2018) or RT-TDDFTRepisky  _et al._ (2015);
Konecny  _et al._ (2018); De Santis  _et al._ (2020); Goings  _et al._ (2016);
Konecny  _et al._ (2016); Kasper  _et al._ (2018, 2020); Darapaneni  _et al._
(2020) should be invoked to account for relativistic effects, which are
sizable already for core excitations even of very light elements. As for
nonrelativistic/relativistic LR-TDDFT and RT-TDDFT themselves, LR-TDDFT is
certainly advantageous over RT-TDDFT in accessing dark states (which are often
important for, e.g., excited energy/charge transfers) as well as in the
assignment of excited states due to the availability of CIS (configuration
interaction singles) type of state eigenvectorsCasida (1995), especially when
full molecular symmetry is employedPeng  _et al._ (2009). Moreover, there
exist very robust algorithmsHuang  _et al._ (2017); Huang and Liu (2019);
Neville and Schuurman (2019) that can access directly the core excited states
as interior roots of the LR-TDDFT eigenvalue problem, without the need to
invoke the so-called core-valence separationStener  _et al._ (2003); Ray  _et
al._ (2007); Besley and Asmuruf (2010); Liang  _et al._ (2011); Van Kuiken
_et al._ (2012), a very good approximation thoughHerbst and Fransson (2020).
Nevertheless, RT-TDDFT has the following advantages over LR-TDDFT: (1) only
the exchange-correlation (XC) potential is needed, requiring no response of
the XC potential, (2) the entire absorption spectrum can be obtained at once,
via a single spectral analysis of the time-dependent dipole moment, (3) the
memory footprint is only about one/two times that of the
relativistic/nonrelativistic ground-state DFT calculation, thereby avoiding
the memory bottleneck of LR-TDDFT resulting from the use of a large number of
trial vectors during the iterative partial diagonalization, and (4) RT-TDDFT
has the capability of accessing dynamical properties beyond the linear
response regime. Because of these, RT-TDDFT has gained great popularity in the
last decades.

From the computational point of view, the very first issue of RT-TDDFT is how
to design a suitable propagator for the numerical integration of the
underlying time-dependent Kohn-Sham equation (TDKS). While a low-order
propagatorKosloff (1988); Castro  _et al._ (2004); Gomez Pueyo  _et al._
(2018) (e.g., the second-order Magus propagator (MP2)Magnus (1954)) is already
sufficient for the low-lying valence absorption spectra, it is no longer the
case for the XAS of systems composed even only of light elements: a low-order
propagator requires an exceedingly small time step to achieve a sufficient
accuracy. The simulation is then much less efficient than the use of a higher-
order propagator along with a larger time step. It is found here that, for the
same spectral accuracy ($<0.05$ eV), the fourth-order optimized commutator-
free exponential time-propagator (oCFET4)Alvermann and Fehske (2011) is most
efficient among the tested high-order propagators (including MPnMagnus (1954);
Blanes  _et al._ (2009), CFETnBlanes and Moan (2006), and oCFETnAlvermann and
Fehske (2011) with n=4, 6) and can hence be regarded as a ‘just good enough’
high-order time propagator for simulations of XAS. However, this is only true
when it is combined with the EPEP3 (exponential prediction of density matrix
and exponential correction of density matrix at three Gauss-Legendre points)
type of predictor-corrector proposed here. Moreover, very robust linear
relations and block aliasing conditions (BAC) for automatic determination of
the largest possible time step as well as a good estimate of the shortest
possible total simulation time will be established to automate the
EPEP3@oCFET4/RT-TDDFT simulations of XAS. Because of the automated nature, the
algorithms will be dubbed collectively as AutoPST (automated determination of
propagator, step, and time), which can further be combined with full molecular
symmetryPeng  _et al._ (2009). To make the presentation complete, we first
recapitulate the TDKS in Sec. II.1 and then discuss the implementation of
oCFET4 in Sec. II.2 and EPEP3 in Sec. II.3, where the closely related
CFET4Blanes and Moan (2006) and EPEP2 are also presented for comparison. The
automated determinations of time step and total time are detailed in Secs.
II.4 and II.5, respectively. Illustrative examples are then provided in Sec.
III to reveal the performance of the proposed AutoPST approach. Some
concluding remarks are finally made in Sec. IV.

##  II Integration of TDKS

###  II.1 TDKS

RT-TDDFT amounts to integrating numerically the TDKS equation, followed by a
spectral analysis to extract information on the electronic transitions. Under
the adiabatic approximation, the TDKS equation discretized in a time-
independent, orthonormal basis $\\{\chi_{\mu}\\}$ can be written as

$i\frac{d}{dt}\mathbf{P}(t)=[\mathbf{F}(t),\mathbf{P}(t)],\quad\mathbf{F}(t)=\mathbf{F}_{0}(t)+\mathbf{V}^{(1)}(t),$
(1)

which is also called Liouville-von Neumann equation for the evolution of the
one-particle density matrix (1PDM) $\mathbf{P}(t)$. Here, $\mathbf{F}_{0}(t)$
is the matrix representation of the field-free KS Hamiltonian at time $t$,

$\displaystyle F_{0}(t)$
$\displaystyle=T+V_{\mathrm{nuc}}+V_{h}[\rho(t)]+c_{\mathrm{x}}\Sigma_{\mathrm{x}}[\gamma(t)]+c_{\mathrm{xc}}V_{\mathrm{xc}}[\rho(t)],$
(2) $\displaystyle\rho(\bm{r},t)$
$\displaystyle=\sum_{\mu\nu}\chi_{\mu}^{\dagger}(\bm{r})\chi_{\nu}(\bm{r})P_{\nu\mu}(t),$
$\displaystyle\gamma(\bm{r},\bm{r}^{\prime},t)$
$\displaystyle=\sum_{\mu\nu}\chi_{\nu}(\bm{r})\chi_{\mu}^{\dagger}(\bm{r}^{\prime})P_{\nu\mu}(t),$
(3)

which is composed of the kinetic energy operator $T$ (nonrelativistic or
relativisticLiu (2020a, b)), (static) nuclear attraction $V_{\mathrm{nuc}}$,
Hartree potential $V_{h}$, exact exchange potential $\Sigma_{\mathrm{x}}$, and
XC potential $V_{\mathrm{xc}}$. The scaling coefficients $c_{\mathrm{x}}$ and
$c_{\mathrm{xc}}$ denote the respective portions of $\Sigma_{\mathrm{x}}$ and
$V_{\mathrm{xc}}$ included in the chosen hybrid functional. As for the time-
dependent external field, we consider here simply a delta-type impulse

$\bm{\varepsilon}(\bm{r},t)=\kappa\bm{n}\delta(t),$ (4)

where $\kappa$ and $\bm{n}$ represent the field strength and polarization unit
vector, respectively. The system-light interaction operator
$\mathbf{V}^{(1)}(t)$ in Eq. (1) then reads (under the dipole approximation)

$\mathbf{V}^{(1)}(t)=\tilde{\mathbf{V}}^{(1)}\delta(t),\quad\tilde{V}^{(1)}_{\mu\nu}=\kappa\sum_{\alpha=x,y,z}\langle\chi_{\mu}|r_{\alpha}|\chi_{\nu}\rangle
n_{\alpha}.$ (5)

Given an initial value $\mathbf{P}(t_{0})$, the solution of Eq. (1) can
formally be written as

$\mathbf{P}(t)=\mathbf{U}(t,t_{0})\mathbf{P}(t_{0})\mathbf{U}^{\dagger}(t,t_{0}),$
(6)

where $\mathbf{U}(t,t_{0})$ takes the following form

$\displaystyle\mathbf{U}(t,t_{0})$
$\displaystyle=\mathcal{T}\mathrm{exp}\\{-i\int_{t_{0}}^{t}\mathbf{F}[\tau]d\tau\\}$
(7)
$\displaystyle=\sum_{n=0}^{\infty}\frac{(-i)^{n}}{n!}\int_{t_{0}}^{t}d\tau_{1}\int_{t_{0}}^{t}d\tau_{2}\cdots\int_{t_{0}}^{t}d\tau_{n}\mathcal{T}[\mathbf{F}(\tau_{1})\mathbf{F}(\tau_{2})\cdots\mathbf{F}(\tau_{n})].$
(8)

Here, $\mathcal{T}$ is the time-ordering operator accounting for the fact that
the KS matrices at different times do not commute. In practice, the time
interval $t-t_{0}$ will be split into $N$ slices with equal length $\Delta t$,
which amounts to splitting $\mathbf{U}(t,t_{0})$ as

$\mathbf{U}(t,t_{0})=\prod^{N-1}_{k=0}\mathbf{U}(t_{k+1},t_{k}),\quad
t_{k}=t_{0}+k\Delta t.$ (9)

In particular, $\mathbf{U}(0^{+},0^{-})$ takes a very simple for the delta-
perturbation (5),

$\mathbf{U}(0^{+},0^{-})=\mathrm{exp}(-i\tilde{\mathbf{V}}^{(1)}),$ (10)

which can readily be obtained from Eq. (7) by noticing that
$\mathbf{F}_{0}(\tau)$ does not change during the infinitely small time
interval $[0^{-},0^{+}]$, such that only the $\delta(t)$ term contributes to
the time integral (see Ref. 57 for a more rigorous derivation). This suggests
that a perturb-then-propagate procedureBaer and Neuhauser (2004); Wang  _et
al._ (2007); Repisky  _et al._ (2015) can be invoked for the evolution of
$\mathbf{P}(t)$, viz.,

$\displaystyle\mathbf{P}(t)$
$\displaystyle=\mathbf{U}(t,0^{+})\mathbf{P}(0^{+})\mathbf{U}^{\dagger}(t,0^{+}),$
(11) $\displaystyle\mathbf{P}(0^{+})$
$\displaystyle=\mathrm{exp}(-i\tilde{\mathbf{V}}^{(1)})\mathbf{P}(0^{-})\mathrm{exp}(i\tilde{\mathbf{V}}^{(1)}),$
(12)

where $\mathbf{P}(0^{-})$ refers to the 1PDM in the absence of the external
field.

Now the crucial point is how to parameterize the local evolution operators
$\mathbf{U}(t_{k+1},t_{k})$. It should be clear from the outset that the
various low-order schemesKosloff (1988); Castro  _et al._ (2004); Gomez Pueyo
_et al._ (2018) are only appropriate for the lowest part of the TDKS spectra,
for they would require exceedingly small time steps to be stable and accurate
enough for the high-end part of the TDKS spectra. After extensive
experimentations, it is learned here that the oCFET4 with three
exponentialsAlvermann and Fehske (2011) performs much better than MP4Magnus
(1954); Blanes  _et al._ (2009) in both accuracy and efficiency. Even higher-
order propagators (MP6Blanes  _et al._ (2009), CFET6Blanes and Moan (2006),
and oCFET6Alvermann and Fehske (2011), etc) are too expensive to be practical.
For comparison, the closely related CFET4 with two exponentialsBlanes and Moan
(2006) is also presented here. Both oCFET4 and CFET4 are unitary by
construction, thereby conserving time reversal symmetry (in the absence of
magnetic fields) as well as the trace, hermiticity, and idempotency of
$\mathbf{P}(t)$. However, strictly speaking, they are appropriate only for
cases where the Hamiltonians governing the dynamics are time independent. For
a time-dependent Hamiltonian as in RT-TDDFT, they have to be combined with
some predictor-correctors to account for the inherent non-autonomy (see Sec.
II.3).

###  II.2 CFET4 and oCFET4 Integrators

The Magnus expansionMagnus (1954); Blanes  _et al._ (2009) has to be
introduced first to make the presentation of CFETBlanes and Moan (2006) clear.
Briefly, the Magnus expansion amounts to replacing the complicated time-
ordered integration

$\displaystyle\mathbf{U}(t_{k+1},t_{k})=\mathcal{T}\mathrm{exp}[\int_{t_{k}}^{t_{k+1}}\mathbf{A}(\tau)d\tau],\quad\mathbf{A}(\tau)=-i\mathbf{F}(\tau),$
(13)

with a time-unordered exponential of an infinite series

$\displaystyle\mathbf{U}(t_{k+1},t_{k})$
$\displaystyle=\mathrm{exp}[\bm{\Omega}(t_{k+1},t_{k})]$ (14)
$\displaystyle=\mathrm{exp}[\sum_{l=1}^{\infty}\varepsilon^{l}\bm{\Omega}_{l}(t_{k+1},t_{k})]_{\varepsilon=1},$
(15)

where the first two wave operators can readily be obtained as

$\displaystyle\bm{\Omega}_{1}(t_{k+1},t_{k})$
$\displaystyle=\int_{t_{k}}^{t_{k+1}}d\tau\mathbf{A}(\tau),$ (16)
$\displaystyle\bm{\Omega}_{2}(t_{k+1},t_{k})$
$\displaystyle=\frac{1}{2}\int_{t_{k}}^{t_{k+1}}d\tau_{1}\int_{t_{k}}^{\tau_{1}}d\tau_{2}[\mathbf{A}(\tau_{1}),\mathbf{A}(\tau_{2})].$
(17)

To evaluate the above time integrals, a Taylor expansion of $\mathbf{A}(\tau)$
can first be made around the midpoint $t_{1/2}=t_{k}+\Delta t/2$,

$\mathbf{A}(\tau)=\sum_{j=0}^{\infty}\mathbf{a}_{j}(\tau-
t_{1/2})^{j},\quad\mathbf{a}_{j}=\frac{1}{j!}\left.\frac{d^{j}\mathbf{A}(\tau)}{d\tau^{j}}\right|_{\tau=t_{1/2}}.$
(18)

Further denoting $\bm{\alpha}_{i}=\mathbf{a}_{i-1}\Delta t^{i}$, it can be
shownBlanes  _et al._ (2009) that only $\\{\bm{\alpha}_{i}\\}_{i=1}^{s}$ are
needed to construct a wave operator $\bm{\Omega}^{[2s]}$ that is correct up to
order $\mathcal{O}(\Delta t^{2s})$. Specifically,

$\displaystyle\bm{\Omega}^{[2]}(t_{k+1},t_{k})$
$\displaystyle=\bm{\alpha}_{1},$ (19)
$\displaystyle\bm{\Omega}^{[4]}(t_{k+1},t_{k})$
$\displaystyle=\bm{\alpha}_{1}-\frac{1}{12}[\bm{\alpha}_{1},\bm{\alpha}_{2}].$
(20)

To facilitate the numerical evaluation of $\\{\bm{\alpha}_{i}\\}_{i=1}^{s}$,
one can first define the following moments of $\mathbf{A}(\tau)$

$\displaystyle\mathbf{A}^{(i-1)}$ $\displaystyle=\frac{1}{\Delta
t^{i-1}}\int_{t_{k}}^{t_{k+1}}(\tau-t_{1/2})^{i-1}\mathbf{A}(\tau)d\tau,\quad
i\in[1,s],$ (21) $\displaystyle=\Delta
t\int_{0}^{1}(\tau-\frac{1}{2})^{i-1}\mathbf{A}(t_{k}+\tau\Delta t)d\tau,$
(22) $\displaystyle=\Delta
t\sum_{j=1}^{n}Q_{ij}^{(s,n)}\mathbf{A}_{j}^{[2s]}+\mathcal{O}(\Delta
t^{2s+1}),$ (23) $\displaystyle\mathbf{A}_{j}^{[2s]}$
$\displaystyle=\mathbf{A}(t_{k}+c_{j}^{[2s]}\Delta t),\quad
Q_{ij}^{(s,n)}=w_{j}^{[2s]}(c_{j}^{[2s]}-\frac{1}{2})^{i-1},$ (24)

where $w_{j}^{[2s]}$ and $c_{j}^{[2s]}$ are the weights and nodes of a
quadrature required to evaluate $\mathbf{A}^{(0)}$ correctly up to order $2s$,
which are applied to the other moments $\mathbf{A}^{(i)}$ as well. Inserting
Eq. (18) into Eq. (21) leads to

$\displaystyle\mathbf{A}^{(i-1)}$
$\displaystyle=\sum_{j=1}^{s}(T^{(s)})_{ij}\bm{\alpha}_{j},$
$\displaystyle(T^{(s)})_{ij}$
$\displaystyle=\frac{1+(-1)^{i+j}}{(i+j-1)2^{i+j-1}},\quad i\in[1,s].$ (25)

Combining Eqs. (23) and (25) then gives rise to the general expression for
$\\{\bm{\alpha}_{j}\\}_{j=1}^{s}$

$\bm{\alpha}_{j}=\Delta t\sum_{k=1}^{n}g_{jk}\mathbf{A}_{k}^{[2s]},\quad
g_{jk}=[(\mathbf{T}^{(s)})^{-1}\mathbf{Q}^{(s,n)}]_{jk},\quad j\in[1,s].$ (26)

Consider the Gauss-Legendre quadrature. It is obvious that only one grid point
is needed for the second-order wave operator $\bm{\Omega}^{[2]}$, i.e.,
$w_{1}^{[2]}=1$ and $c_{1}^{[2]}=\frac{1}{2}$. We then have $Q^{(1,1)}=1$ and
$(T^{(1)})^{-1}=1$ and hence

$\bm{\Omega}^{[2]}(t_{k+1},t_{k})=\bm{\alpha}_{1}=\mathbf{A}(t_{k}+\frac{1}{2}\Delta
t)\Delta t,$ (27)

which gives rise to the standard, second-order midpoint Magnus propagator
(MP2)

$\mathbf{U}^{[2]}(t_{k+1},t_{k})=\mathrm{exp}[\bm{\Omega}^{[2]}(t_{k+1},t_{k})].$
(28)

In contrast, the evaluation of $\bm{\Omega}^{[4]}$ requires two grid points,
i.e., $w_{1}^{[4]}=\frac{1}{2}$, $c_{1}^{[4]}=\frac{1}{2}-\frac{\sqrt{3}}{6}$,
$w_{2}^{[4]}=\frac{1}{2}$, and $c_{2}^{[4]}=\frac{1}{2}+\frac{\sqrt{3}}{6}$,
which give rise to

$\mathbf{Q}^{(2,2)}=\begin{pmatrix}\frac{1}{2}&\frac{1}{2}\\\
-\frac{\sqrt{3}}{12}&\frac{\sqrt{3}}{12}\end{pmatrix},\quad(\mathbf{T}^{(2)})^{-1}=\begin{pmatrix}1&0\\\
0&12\end{pmatrix},$ (29)

and hence

$\displaystyle\begin{pmatrix}\bm{\alpha}_{1}\\\ \bm{\alpha}_{2}\end{pmatrix}$
$\displaystyle=\Delta
t(\mathbf{T}^{(2)})^{-1}\mathbf{Q}^{(2,2)}\begin{pmatrix}\mathbf{A}_{1}^{[4]}\\\
\mathbf{A}_{2}^{[4]}\end{pmatrix}$ $\displaystyle=\Delta
t\begin{pmatrix}\frac{1}{2}(\mathbf{A}_{1}^{[4]}+\mathbf{A}_{2}^{[4]})\\\
\sqrt{3}(\mathbf{A}_{2}^{[4]}-\mathbf{A}_{1}^{[4]})\end{pmatrix},$ (30)
$\displaystyle\bm{\Omega}^{[4]}(t_{k+1},t_{k})$
$\displaystyle=\frac{1}{2}(\mathbf{A}_{1}^{[4]}+\mathbf{A}_{2}^{[4]})\Delta
t-\frac{\sqrt{3}}{12}[\mathbf{A}_{1}^{[4]},\mathbf{A}_{2}^{[4]}]\Delta t^{2}.$
(31)

The fourth-order Magnus propagator (MP4) then reads

$\displaystyle\mathbf{U}^{[4]}(t_{k+1},t_{k})$
$\displaystyle=\mathrm{exp}[\bm{\Omega}^{[4]}(t_{k+1},t_{k})].$ (32)

Instead of the above Magnus expansion which involves increasingly nested
commutators at high orders, the commutator-free type of expansion of Eq. (13)
assumes the following form

$\displaystyle\mathbf{U}^{[2s]}_{\mathrm{CF},m}(t_{k+1},t_{k})=e^{\widetilde{\bm{\Omega}}}=e^{\widetilde{\bm{\Omega}}_{m}}e^{\widetilde{\bm{\Omega}}_{m-1}}\cdots
e^{\widetilde{\bm{\Omega}}_{1}},$ (33)
$\displaystyle\widetilde{\bm{\Omega}}_{i}=\sum_{j=1}^{s}f_{ij}\bm{\alpha}_{j}=\Delta
t\sum_{k=1}^{n}x_{ik}\mathbf{A}_{k}^{[2s]},\quad
x_{ik}=\sum_{j=1}^{s}f_{ij}g_{jk},$ (34)

where $m$ denotes the number of exponentials. The expansion coefficients
$f_{ij}$ can be determinedBlanes and Moan (2006) by equating the terms of
$\widetilde{\bm{\Omega}}$ to those of the Magnus expansion $\bm{\Omega}$ (14)
at each order, up to $\mathcal{O}(\Delta t^{2s})$. Consider CFET4 with two
exponentials

$\mathbf{U}_{\mathrm{CF},2}^{[4]}(t_{k+1},t_{k})=e^{f_{21}\bm{\alpha}_{1}+f_{22}\bm{\alpha}_{2}}e^{f_{11}\bm{\alpha}_{1}+f_{12}\bm{\alpha}_{2}},$
(35)

where the coefficients should be subject to

$f_{21}=f_{11},\quad f_{22}=-f_{12}$ (36)

to fulfill time reversal symmetry. Further use of the Baker-Campbell-Hausdorff
formula can be made to merge the two exponentials into one

$\mathbf{U}_{\mathrm{CF},2}^{[4]}(t_{k+1},t_{k})=e^{2f_{11}\bm{\alpha}_{1}+\frac{f_{12}}{2}[\bm{\alpha}_{1},\bm{\alpha}_{2}]},$
(37)

the exponent of which can be equated to $\bm{\Omega}^{[4]}(t_{k+1},t_{k})$ in
Eq. (20), thereby leading to

$\mathbf{U}_{\mathrm{CF},2}^{[4]}(t_{k+1},t_{k})=e^{\frac{1}{2}\bm{\alpha}_{1}+\frac{1}{6}\bm{\alpha}_{2}}e^{\frac{1}{2}\bm{\alpha}_{1}-\frac{1}{6}\bm{\alpha}_{2}},$
(38)

the proper implementation of which is described in Algorithm 1 in Appendix A.

Instead of the termwise matching between CFETn and MPn, the expansion
coefficients $\\{f_{ij}\\}$ in Eq. (34) can also be optimizedAlvermann and
Fehske (2011) to minimize the error between the approximate
$\widetilde{\bm{\Omega}}$ and the exact $\bm{\Omega}$, thereby leading to the
optimized (oCFET) variants of CFET. Without going into further details, the
required coefficients $\\{x_{ik}\\}$ can be determined by reexpressing the
coefficients $\mathbf{A}_{n}$ defined in equation (25) of Ref. 75 in terms of
the present moments $\mathbf{A}_{j}^{[2s]}$ (24). The oCFET4 with three
exponentials considered here is specified in Algorithm 2 in Appendix A.

###  II.3 Predictor-Correctors

The above CFET4 and oCFET4 have to be combined with suitable predictor-
correctors to account for the time dependence in RT-TDDFT. The available
predictor-correctors Castro  _et al._ (2004); Cheng  _et al._ (2006); Zhu and
Herbert (2018); Repisky  _et al._ (2015) were designed only for MP2 and hence
cannot be used here. Instead, we got to design new predictor-correctors for
CFET4 and oCFET4 case by case.

####  II.3.1 Predictor-Corrector for CFET4

As a fourth-order propagator, CFET4 requires two Gauss-Legendre quadrature
points. For this reason, an EPEP2 (exponential prediction of density matrix
and exponential correction of density matrix at two Gauss-Legendre points)
type of predictor-corrector is proposed here, see Fig. 1 (and Algorithm 3 in
Appendix A). The workflow goes as follows.

Step (1): Propagate $\mathbf{P}(t_{k})$ to $\mathbf{P}(t_{1}^{[4]})$ with the
second-order enforced time reversal symmetry (ETRS) propagatorCastro  _et al._
(2004) (see Algorithm 4 in Appendix A),
$\displaystyle\mathbf{U}_{\mathrm{ETRS}}(t_{k+1},t_{k})$
$\displaystyle=e^{-i\frac{\Delta t}{2}\mathbf{F}(t_{k+1})}e^{-i\frac{\Delta
t}{2}\mathbf{F}(t_{k})},$ (39) $\displaystyle\mathbf{P}(t_{1}^{[4]})$
$\displaystyle=\mathbf{U}_{\mathrm{ETRS}}(t_{k+1},t_{k})\mathbf{P}(t_{k})\mathbf{U}_{\mathrm{ETRS}}^{\dagger}(t_{k+1},t_{k}),$
$\displaystyle t_{1}^{[4]}$ $\displaystyle=t_{k}+c_{1}^{[4]}\Delta t.$ (40)
Step (2): Construct $\mathbf{F}(t_{1}^{[4]})$ with $\mathbf{P}(t_{1}^{[4]})$.
Step (3): Propagate $\mathbf{P}(t_{1}^{[4]})$ to $\mathbf{P}(t_{2}^{[4]})$
with Euler’s rule,  $\displaystyle\mathbf{P}(t_{2}^{[4]})$
$\displaystyle=e^{-i\mathbf{F}(t_{1}^{[4]})(t_{2}^{[4]}-t_{1}^{[4]})}\mathbf{P}(t_{1}^{[4]})e^{i\mathbf{F}(t_{1}^{[4]})(t_{2}^{[4]}-t_{1}^{[4]})},$
$\displaystyle t_{2}^{[4]}$ $\displaystyle=t_{k}+c_{2}^{[4]}\Delta t.$ (41)
Step (4): Construct $\mathbf{F}(t_{2}^{[4]})$ with $\mathbf{P}(t_{2}^{[4]})$.
Step (5): Propagate $\mathbf{P}(t_{k})$ to $\mathbf{P}(t_{k+1})$ with CFET4
(see Algorithm 1 in Appendix A).  Step (6): Construct $\mathbf{F}(t_{k+1})$
with $\mathbf{P}(t_{k+1})$.  Step (7): Backward propagate
$\mathbf{P}(t_{k+1})$ to $\mathbf{P}(t_{2}^{[4]})$ with Euler’s rule,
$\mathbf{P}(t_{2}^{[4]})=e^{i\mathbf{F}(t_{k+1})(t_{k+1}-t_{2}^{[4]})}\mathbf{P}(t_{k+1})e^{-i\mathbf{F}(t_{k+1})(t_{k+1}-t_{2}^{[4]})},$
(42) so as to correct $\mathbf{P}(t_{2}^{[4]})$.  Step (8): Repeat steps (4)
to (7) until $\mathbf{P}(t_{k+1})$ converges according to the following
criterionZhu and Herbert (2018)
$\frac{\|\mathbf{X}_{1}-\mathbf{X}_{2}\|_{F}}{M}<\xi,$ (43) where
$\|\cdot\|_{F}$ denotes the Frobenius norm, $\mathbf{X}_{1}$ and
$\mathbf{X}_{2}$ represents two $M\times M$ matrices to be compared, and $\xi$
is the chosen threshold.

In the above procedure, $\mathbf{P}(t_{1}^{[4]})$ does not enter the iteration
cycle (steps (4) to (7)) but is predicted directly from $\mathbf{P}(t_{k})$.
That is, $\mathbf{P}(t_{1}^{[4]})$ is not to be corrected iteratively by
$\mathbf{P}(t_{k+1})$ although this can in principle be done. Such an option
arises from the following consideration: $t_{1}^{[4]}$ is close to $t_{k}$,
such that the second-order ETRS (instead of the first-order Euler) propagator
would ensure an accurate prediction of $\mathbf{P}(t_{1}^{[4]})$. In contrast,
had $\mathbf{P}(t_{1}^{[4]})$ been corrected by $\mathbf{P}(t_{k+1})$, the
accuracy would be deteriorated due to the larger interval between
$t_{1}^{[4]}$ and $t_{k+1}$.

Figure 1: One time-step flowchart of EPEP2 for CFET4. In circles are step
indices of Algorithm 3 in Appendix A.

####  II.3.2 Predictor-Corrector for oCFET4

Albeit a fourth-order propagator, oCFET4 requires three Gauss-Legendre
quadrature points. In the same spirit as EPEP2, an EPEP3 type of predictor-
corrector is proposed here, see Fig. 2 (and Algorithm 5 in Appendix A). The
workflow goes as follows.

Step 1: Propagate $\mathbf{P}(t_{k})$ to $\mathbf{P}(t_{1}^{[6]})$ with ETRS.
Here, $t_{1}^{[6]}=t_{k}+c_{1}^{[6]}\Delta t$ and
$c_{1}^{[6]}=\frac{1}{2}-\frac{\sqrt{15}}{10}$.  Step 2: Construct
$\mathbf{F}(t_{1}^{[6]})$ with $\mathbf{P}(t_{1}^{[6]})$.  Step 3: Propagate
$\mathbf{P}(t_{1}^{[6]})$ to $\mathbf{P}(t_{2}^{[6]})$ with Euler’s rule,
$\displaystyle\mathbf{P}(t_{2}^{[6]})$
$\displaystyle=e^{-i\mathbf{F}(t_{1}^{[6]})(t_{2}^{[6]}-t_{1}^{[6]})}\mathbf{P}(t_{1}^{[6]})e^{i\mathbf{F}(t_{1}^{[6]})(t_{2}^{[6]}-t_{1}^{[6]})},$
$\displaystyle t_{2}^{[6]}$ $\displaystyle=t_{k}+c_{2}^{[6]}\Delta t,\quad
c_{2}^{[6]}=\frac{1}{2}.$ (44) Step 4: Construct $\mathbf{F}(t_{2}^{[6]})$
with $\mathbf{P}(t_{2}^{[6]})$.  Step 5: Propagate $\mathbf{P}(t_{1}^{[6]})$
to $\mathbf{P}(t_{3}^{[6]})$ with MP2 (28),
$\displaystyle\mathbf{P}(t_{3}^{[6]})$
$\displaystyle=e^{-i\mathbf{F}(t_{2}^{[6]})(t_{3}^{[6]}-t_{1}^{[6]})}\mathbf{P}(t_{3}^{[6]})e^{i\mathbf{F}(t_{2}^{[6]})(t_{3}^{[6]}-t_{1}^{[6]})},$
$\displaystyle t_{3}^{[6]}$ $\displaystyle=t_{k}+c_{3}^{[6]}\Delta t,\quad
c_{3}^{[6]}=\frac{1}{2}+\frac{\sqrt{15}}{10}.$ (45) Step 6: Construct
$\mathbf{F}(t_{3}^{[6]})$ with $\mathbf{P}(t_{3}^{[6]})$.  Step 7: Propagate
$\mathbf{P}(t_{k})$ to $\mathbf{P}(t_{k+1})$ with oCFET4 (see Algorithm 2 in
Appendix A).  Step 8: Construct $\mathbf{F}(t_{k+1})$ with
$\mathbf{P}(t_{k+1})$.  Step 9: Backward propagate $\mathbf{P}(t_{k+1})$ to
$\mathbf{P}(t_{3}^{[6]})$ with Euler’s rule,
$\mathbf{P}(t_{3}^{[6]})=e^{i\mathbf{F}(t_{k+1})(t_{k+1}-t_{3}^{[6]})}\mathbf{P}(t_{k+1})e^{-i\mathbf{F}(t_{k+1})(t_{k+1}-t_{3}^{[6]})},$
(46) so as to correct $\mathbf{P}(t_{3}^{[6]})$.  Step 10: Construct
$\mathbf{F}(t_{3}^{[6]})$ with $\mathbf{P}(t_{3}^{[6]})$.  Step 11: Backward
propagate $\mathbf{P}(t_{3}^{[6]})$ to $\mathbf{P}(t_{4}^{[6]})$ with Euler’s
rule,  $\displaystyle\mathbf{P}(t_{4}^{[6]})$
$\displaystyle=e^{i\mathbf{F}(t_{3}^{[6]})(t_{3}^{[6]}-t_{4}^{[6]})}\mathbf{P}(t_{3}^{[6]})e^{-i\mathbf{F}(t_{3}^{[6]})(t_{3}^{[6]}-t_{4}^{[6]})},$
$\displaystyle t_{4}^{[6]}$ $\displaystyle=\frac{t_{2}^{[6]}+t_{3}^{[6]}}{2}.$
(47) Step 12: Construct $\mathbf{F}(t_{4}^{[6]})$ with
$\mathbf{P}(t_{4}^{[6]})$.  Step 13: Backward propagate
$\mathbf{P}(t_{3}^{[6]})$ to $\mathbf{P}(t_{2}^{[6]})$ with MP2 (28),
$\mathbf{P}(t_{2}^{[6]})=e^{i\mathbf{F}(t_{4})(t_{3}^{[6]}-t_{2}^{[6]})}\mathbf{P}(t_{3}^{[6]})e^{-i\mathbf{F}(t_{4})(t_{3}^{[6]}-t_{2}^{[6]})},$
(48) so as to correct $\mathbf{P}(t_{2}^{[6]})$.  Step 14: Construct
$\mathbf{F}(t_{2}^{[6]})$ with $\mathbf{P}(t_{2}^{[6]})$.  Step 15: Repeat
steps (7) to (14) until $\mathbf{P}(t_{k+1})$ converges according to criterion
(43).

Figure 2: One time-step flowchart of EPEP3 for oCFET4. In circles are step
indices of Algorithm 5 in Appendix A.

In the above procedure, $\mathbf{P}(t_{1}^{[6]})$ is predicted directly from
$\mathbf{P}(t_{k})$, for $t_{1}^{[6]}$ is much closer to $t_{k}$ than to
$t_{k+1}$. This stays in the same spirit as EPEP2. $\mathbf{P}(t_{3}^{[6]})$
is corrected by using the backward Euler propagation of $\mathbf{P}(t_{k+1})$,
for $t_{3}^{[6]}$ is closer to $t_{k+1}$. Since the interval
$t_{3}^{[6]}-t_{2}^{[6]}$ for $\mathbf{P}(t_{2}^{[6]})$ is larger than the
interval $t_{k+1}-t_{3}^{[6]}$ for $\mathbf{P}(t_{3}^{[6]})$, an auxiliary
time moment $t_{4}^{[6]}$ has been introduced to correct
$\mathbf{P}(t_{2}^{[6]})$ by using the backward MP2 propagation from
$\mathbf{P}(t_{3}^{[6]})$.

Since EPEP3 and EPEP2 share the same strategy, they should converge at the
same rate. As a matter of fact, in the present calculations of core
excitations with the predicted time steps (see Sec. II.4), the first iteration
already gives rise to accurate predictions of $\mathbf{P}(t_{k+1})$, such that
the second iteration is just to confirm the convergence, which is far below
the chosen threshold $\xi=10^{-7}$ in Eq. (43).

###  II.4 Time Step

Having chosen an appropriate propagator (e.g., MP2 for low-lying valence
excitations and oCFET4 for core excitations), along with the corresponding
corrector-predictor, the next question is how to choose an appropriate time
step $\Delta t$. Instead of playing around with $\Delta t$ in a trial-and-
error manner, we here try to fix it in an _a priori_ manner, so as to
facilitate the dynamics simulations. It is well known that $\Delta t$ is
bounded by both the stability and accuracy of the chosen propagator as well as
the Nyquist criterion for signal samplingNyquist (1928); Shannon (1948), i.e.,
$\Delta t\leq\Delta t_{H}=\frac{1}{f_{s}}<\frac{1}{2f_{H}}$, with $f_{H}$
being the highest temporal frequency contained in the sampled signal.
Literally, the sampling frequency $f_{s}$ ($=1/\Delta t$) must be twice larger
than the width $f_{H}-(-f_{H})$ of the spectrum contained in the band limited
real signal. This restriction is to avoid the so-called aliasing phenomenon,
which stems from the fact that every frequency within $[-f_{H},f_{H}]$ will be
copied infinitely many times (separated by the sampling frequency $f_{s}$) by
the Fourier transform of the sampled signal. However, if the spectrum
$[-f_{H},f_{H}]$ is separated in blocks and only one particular block is of
interest, an even larger time step $\Delta t$ can be chosen to sample the
signal. To see this, consider the spectrum of a real-valued time signal shown
in Fig. 3, which has four blocks, $[-f_{H},-f_{L}]$, $[-B,0]$, $[0,B]$, and
$[f_{L},f_{H}]$. To extract the frequencies within $[0,B]$, the proper
sampling frequency $f_{s}$ can be determined by the requirement that no alias
from the other blocks will appear in the interval $[0,B]$. Consider first the
block $[f_{L},f_{H}]$, which is shifted to the left by $(n-1)$ times until the
right end of $[0,B]$. The next left shift of $[f_{L},f_{H}]$ must cross
$[0,B]$ in one step of $f_{s}$, thereby resulting in the first necessary
condition for $f_{s}$,

$f_{s}>B+f_{H}-f_{L}.$ (49)

The same is obtained by right shifts of $[-f_{H},-f_{L}]$ in steps of $f_{s}$.
Likewise, the one-step right shift of $[-B,0]$ gives rise to the following
necessary condition

$f_{s}>2B.$ (50)

The two conditions (49) and (50) can be combined to yield the lower bound
$\tilde{f}_{s}$ of $f_{s}$

$\tilde{f}_{s}=\max\\{B+f_{H}-f_{L},2B\\}<f_{s}.$ (51)

Moreover, it is obvious that the following conditions

$\displaystyle f_{L}-(n-1)f_{s}>B,\quad f_{H}-\mathrm{nf}_{s}<0,\quad
n=1,2,\cdots,$ (52)

should also be satisfied simultaneously for the left shifts of
$[f_{L},f_{H}]$. The conditions in Eqs. (51) and (52) can be combined to yield

$\displaystyle\frac{f_{H}}{n}<f_{s}<\frac{f_{L}-B}{n-1},\quad 1\leqslant
n\leqslant\left\lfloor\frac{f_{L}-B}{\tilde{f}_{s}}\right\rfloor+1,$ (53)

where the symbol $\lfloor y\rfloor$ denotes the largest integer less than or
equal to $y$. Similarly, the following conditions

$\displaystyle-f_{L}+(m-1)f_{s}<0,\quad-f_{H}+mf_{s}>B,\quad m=1,2,\cdots,$
(54)

should be satisfied additionally for the right shifts of $[-f_{H},-f_{L}]$.
The conditions in Eqs. (51) and (54) can be combined to yield

$\displaystyle\frac{f_{H}+B}{m}<f_{s}<\frac{f_{L}}{m-1},\quad 1\leqslant
m\leqslant\left\lfloor\frac{f_{L}}{\tilde{f}_{s}}\right\rfloor+1.$ (55)

Finally, no additional condition is required for the right shift of $[-B,0]$.

In summary, to read out correctly the frequencies in the interval $[0,B]$ from
the Fourier transform of the sampled signal, the sampling frequency $f_{s}$
has to satisfy the block aliasing conditions (BAC) (51), (53), and (55)
simultaneously. The analysis can readily be generalized to extract the
frequencies of any block of a spectrum composed of $2N$ blocks. The algorithm
goes as follows.

  1. (1)

For each symmetry $\Gamma$ of the excited states, determine the corresponding
symmetries of the virtual ($\gamma_{p}$) and occupied ($\gamma_{h}$) canonical
molecular orbitals (CMOs) and calculate their energy differences
$\\{\Delta\epsilon_{ai}\\}|_{a\in\gamma_{p},i\in\gamma_{h}}$. Sort
$\\{\Delta\epsilon_{ai}\\}$ of all the $\gamma_{p}$-$\gamma_{h}$ pairs in
ascending order, leading to
$\\{\Delta\epsilon_{n}\\}_{n=1}^{\Gamma,N_{\Gamma}}$.

  2. (2)

Partition $\\{\Delta\epsilon_{n}\\}_{n=1}^{\Gamma,N_{\Gamma}}$ into blocks
$\lceil\Delta\epsilon\rceil^{\Gamma,M_{\Gamma}}_{m=1}$ according to their
gaps. Specifically, if the separation between $\Delta\epsilon_{n}$ and
$\Delta\epsilon_{n+1}$ is larger than $\omega_{G}=2\omega_{b}+\omega_{c}$, set
$\Delta\epsilon_{n}$ and $\Delta\epsilon_{n+1}$ as the uppermost and lowermost
elements of the current and next blocks, respectively. Here, $\omega_{b}$ is a
parameter (e.g., 0.5 a.u.) to broaden each block, so as to account for the
uncertainties of the independent particle approximation (IPA) of the true
excitation energies, whereas $\omega_{c}$ is another parameter (e.g., 0.2
a.u.) to ensure the minimal separation between two adjacent broadened blocks
(see Fig. S1/S2 for an example). The in total $M_{\Gamma}$ blocks
$\frac{1}{2\pi}\lceil\Delta\epsilon\rceil^{\Gamma,M_{\Gamma}}_{m=1}\equiv[f_{Lm},f_{Hm}]_{m=1}^{\Gamma,M_{\Gamma}}$
along with their negative counter parts
$[-f_{Hm},-f_{Lm}]_{m=1}^{\Gamma,M_{\Gamma}}$ form the estimated structure of
the whole TDKS spectrum .

  3. (3)

Choose a target block, say $[f_{Lk},f_{Hk}]$, and calculate the minimal
sampling frequency

$\tilde{f}_{s}=f_{Hk}-f_{Lk}+\max_{m\in[1,M_{\Gamma}]}(f_{Hm}-f_{Lm}).$ (56)

  4. (4)

Build up the corresponding inequalities for each block $[f_{a},f_{b}]$
($f_{b}>f_{a}$; including those of negative frequencies). For those blocks
lower and higher than $[f_{Lk},f_{Hk}]$, we have

$\displaystyle\frac{f_{Hk}-f_{a}}{n}<f_{s}<\frac{f_{Lk}-f_{b}}{n-1},\quad
1\leqslant
n\leqslant\left\lfloor\frac{f_{Lk}-f_{b}}{\tilde{f}_{s}}\right\rfloor+1,$ (57)

and

$\displaystyle\frac{f_{b}-f_{Lk}}{n}<f_{s}<\frac{f_{a}-f_{Hk}}{n-1},\quad
1\leqslant
n\leqslant\left\lfloor\frac{f_{a}-f_{Hk}}{\tilde{f}_{s}}\right\rfloor+1,$ (58)

respectively.

  5. (5)

Determine the intersection of the inequalities of all blocks so as to obtain
the allowed sampling frequency $f_{s}$ and hence the time step $\Delta
t=1/f_{s}$ for each symmetry $\Gamma$ of excited states.

Note that when the whole spectrum is regarded as one big block
$[-f_{H},f_{H}]$, the above algorithm reduces naturally to the Nyquist time
step $\Delta t_{H}=\frac{1}{f_{s}}=\frac{1}{2f_{H}}$, as can be seen from Eq.
(56).

Figure 3:  (a) Schematic representation of a sparse spectrum, with frequencies
located in blocks $[-f_{H},-f_{L}]$, $[-B,0]$, $[0,B]$, and $[f_{L},f_{H}]$.
(b) Green blocks with dashed contour lines are periodic images (aliases) of
the original frequencies in $[f_{L},f_{H}]$, separated by the sampling
frequency $f_{s}$. The conditions in Eq. (53) ensure that the $n$-th and
$(n-1)$-th image blocks of $[f_{L},f_{H}]$ are located on the left and right
hands of sides of $[0,B]$, respectively. Further incorporation of the
conditions in Eq. (55) will separate the $n$-th image block of $[f_{L},f_{H}]$
from $[-B,0]$, as a result of the separation between $[0,B]$ and the $n$-th
image block of $[-f_{H},-f_{L}]$ (because $[-f_{H},-f_{L}]$ and
$[f_{L},f_{H}]$ are symmetric with respect to the 0 Hz axis).

At this stage, it must be realized that the time step determined by the above
algorithm only guarantees that the target frequencies are not messed up with
aliases from the other blocks. The accuracy of the spectral analysis depends
further on the quality of the input time signal. For a given propagator, the
smaller the time step, the more accurate the time signal. On the other hand,
for a given time step, the more accurate the propagator, the more accurate the
time signal. However, a very small time step or a very accurate propagator
usually implies a very high computational cost. Therefore, it is important to
balance these factors in the allowed range of sampling frequency. To achieve
this, we first tried using the following parameter

$\tau_{P}=\frac{\|\mathbf{P}(t)-\mathbf{P}_{\mathrm{ref}}(t)\|_{F}}{M\cdot t}$
(59)

to quantify the accuracy of $\mathbf{P}(t)$ (evolved by the chosen propagator
with time step $\Delta t$) against a reference 1PDM
$\mathbf{P}_{\mathrm{ref}}(t)$ (obtained by CF6:5Opt (see Table 6 in Ref.75)
with $\Delta t_{H}/2$). However, preliminary experimentation reveals that, for
core excitations, $\tau_{P}$ is rather sensitive to all the ingredients
(atomic number, XC functional, electron kinematics, and basis set), meaning
that it is not directly connected to the target accuracy. Having realized that
such ingredients have been ‘renormalized’ into the lowest core level
$\epsilon_{1s}$ of the ground-state calculation, we switch to seek for a
relation between $\Delta t$ and $\epsilon_{1s}$ for core excitations. As can
be seen from Fig. 4, for the lowest five bright $K$-edge core excited states
of HX (X = F, Cl, Br), there exist very good linear relations between
$log_{10}(\Delta t)$ and $log_{10}(|\epsilon_{1s}|)$ for both CFET4 and
oCFET4. Here, $\Delta t$ is the largest time step to guarantee the desired
spectral accuracy of less than 0.05 eV. The fitting functions for different XC
functionals are further documented in Table 1. Noticeably, the slopes and
interceptions are slightly different for different XC functionals but which
are hardly surprising, for different functionals leads to different degrees of
noncommutativity between KS matrices at different times: the more exact
exchange, the more enhanced noncommutativity and hence the smaller the time
step. Although such relations have been obtained based on nonrelativistic
calculations with a single basis set (DEF2-SVPDWeigend and Ahlrichs (2005);
Rappoport and Furche (2010)), they should be applicable to other basis sets
and other systems with relativistic effects accounted for. The underlying
reasons are twofold: (1) different basis sets produce essentially the same
energies for all the occupied and low-lying virtual orbitals and (2)
relativity is heavily dominated by the time-independent one-body terms (which
do not cause noncommutativity between KS matrices at different times). In
other words, such linear relations are very robust with respect to atomic
number, electron kinematics, and basis set, and can hence be applied to all
the elements in the Periodic Table, to achieve an accuracy better than 0.05 eV
for core excitation energies as compared to those by LR-TDDFT. To confirm
this, the time steps predicted by the linear relations in Table 1 (which
requires only the ground-state lowest core orbital energy $\epsilon_{1s}$ as
input) will be taken in Sec. III for the $K$\- and $L$-edge core excitations
of HI and PbO. Note in passing that the time step predicted by the linear
relations should further be verified to comply with the above BACs. If not, it
should be reduced to the largest allowed value.

After having determined the time step for each propagator, it is attempting to
know which propagator will be most efficient in terms of the number of Fock
builds. To this end, the numbers ($N_{\mathrm{FBS}}$) of Fock builds per time
step are summarized in Table 2 for MP2, MP4, CFET4, and oCFET4.

---
(a)
(b)
Figure 4: Linear relations between $log_{10}(\Delta t)$ and
$log_{10}(|\epsilon_{1s}|)$ for the lowest five bright $K$-edge core excited
states of HF, HCl, and HBr by (a) CFET4@RT-TDDFT and (b) oCFET4@RT-TDDFT with
different exchange-correlation functionals (nonrelativistic calculations with
the DEF2-SVPD basis setWeigend and Ahlrichs (2005); Rappoport and Furche
(2010)). $\Delta t$: the largest time step for an accuracy of $<0.05$ eV;
$\epsilon_{1s}$: the ground-state lowest core orbital energy. Table 1: Fitting
functions∗ $log_{10}(\Delta t)=a\cdot log_{10}(|\epsilon_{1s}|)+b$ for
$K$-edge core excitations by CFET4@RT-TDDFT and oCFET4@RT-TDDFT (cf. Fig. 4)

propagator | functional | $a$ | $b$
---|---|---|---
CFET4 |  BLYPBecke (1988); Lee  _et al._ (1988) | -0.96018601 | 0.63991532
B3LYPBecke (1993); Stephens  _et al._ (1994) | -0.96592829 | 0.37019600
LC-BLYPIikura  _et al._ (2001) | -0.97307440 | 0.24557525
CAM-B3LYPYanai  _et al._ (2004) | -0.98942182 | 0.30841338
BHandHLYPBecke (1993, 1988); Lee  _et al._ (1988) | -0.97407896 | 0.26693596
HF | -0.92526295 | 0.02526600
oCFET4 |  BLYPBecke (1988); Lee  _et al._ (1988) | -0.96171191 | 0.81185106
B3LYPBecke (1993); Stephens  _et al._ (1994) | -0.96656731 | 0.64711096
LC-BLYPIikura  _et al._ (2001) | -0.99080653 | 0.62700701
CAM-B3LYPYanai  _et al._ (2004) | -0.99061081 | 0.63794424
BHandHLYPBecke (1993, 1988); Lee  _et al._ (1988) | -0.97647005 | 0.59990197
HF | -0.95558724 | 0.48821463

  * *

$\epsilon_{1s}$: the ground-state lowest core orbital energy.

Table 2: Numbers ($N_{\mathrm{FBS}}$) of Fock builds per time step involved in
the predictor-correctors

Predictor-Corrector | Propagator |  $N_{\mathrm{FBS}}$a | Formal order
---|---|---|---
EPEP2b | CFET4, MP4 | $3\cdot N_{\mathrm{it}}$ | 4
EPEP3c | oCFET4 | $5\cdot N_{\mathrm{it}}-1$ | 4

  * a

$N_{\mathrm{it}}$: average number of iterations in each for loop within a time
step.

  * b

Algorithm 3 in Appendix A.

  * c

Algorithm 5 in Appendix A.

###  II.5 Total Sampling Time

The last ingredient to be specified for the AutoPST approach is the total
sampling time $T$, which depends not only on the desired spectral resolution
but also on the chosen spectral analysis method. The most well-known spectral
analysis technique is the (discrete) Fourier transform, which can extract all
the frequencies contained in the signal at once but at the expense of a very
long simulation time ($T\geq 2\pi/\delta\omega$, with $\delta\omega$ being the
smallest frequency separation). When the spectrum is block-structured (which
is usually the case), it is more efficient to employ techniquesBruner  _et
al._ (2016); Wall and Neuhauser (1995); Mandelshtam and Taylor (1997); Pang
_et al._ (1998); Schelter and Kümmel (2018) that can extract a particular
block of the spectrum. Here the filter diagonalization (FD)
approachMandelshtam and Taylor (1997) is utilized for the $K$\- and $L$-edge
core excitations. In essence, FD is a nonlinear fitting process to obtain $2K$
($\leq N$) unknowns $\\{d_{k},\omega_{k}\\}_{k=1}^{K}$ from $N$ data
$\\{C(t_{n})\\}_{n=0}^{N-1}$, such that

$C(t_{n})=\sum_{k=1}^{K}d_{k}e^{-i\omega_{k}t_{n}},\quad n\in[0,N-1].$ (60)

From the informational point of view, the sampling time $T$ required by FD is
at least $4\pi/\delta\bar{\omega}$ (i.e., $2\pi/\delta\bar{\omega}$ for
$\\{d_{k}\\}_{k=1}^{K}$ and $2\pi/\delta\bar{\omega}$ for
$\\{\omega_{k}\\}_{k=1}^{K}$), where $\delta\bar{\omega}$ is the average
frequency separation in the target block of the spectrum (which is usually
much larger than $\delta\omega$). The lower bound holds when the intrinsic
frequencies $\\{\omega_{k}\\}_{k=1}^{K}$ are uniformly separated. In case that
two frequencies are very close to each other, a time longer than
$4\pi/\delta\bar{\omega}$ will be needed even if $\delta\bar{\omega}$ is not
that small. To incorporate such features of FD, we design the following
algorithm for the determination of $T$:

  1. (a)

For each symmetry $\Gamma$ of the excited states, determine the corresponding
symmetries of the virtual ($\gamma_{p}$) and occupied ($\gamma_{h}$) CMOs and
calculate their energy differences
$\\{\Delta\epsilon_{ai}\\}|_{a\in\gamma_{p},i\in\gamma_{h}}$. If
$|\Delta\epsilon_{m}$-$\Delta\epsilon_{n}|<\omega_{d}$ (e.g.,
$\omega_{d}=0.02$ a.u.), the $m$th and $n$th particle-hole pairs are
considered as degenerate, such that only one of them is retained. Sort
$\\{\Delta\epsilon_{ai}\\}$ of all the $\gamma_{p}$-$\gamma_{h}$ pairs in
ascending order, leading to
$\\{\Delta\epsilon_{n}\\}_{n=1}^{\Gamma,N_{\Gamma}}$.

  2. (b)

Partition $\\{\Delta\epsilon_{n}\\}_{n=1}^{\Gamma,N_{\Gamma}}$ into blocks
$\lfloor\Delta\epsilon\rfloor^{\Gamma,H_{\Gamma}}_{I=1}$ according to their
homogeneity. Specifically, when inserting $\Delta\epsilon_{n+1}$ into the
current block renders the standard deviation of adjacent separation
$\sigma_{\delta\omega}$ larger than the homogeneity $\omega_{\sigma}$ (e.g.,
$0.1$ a.u.), set $\Delta\epsilon_{n}$ as the uppermost and lowermost elements
of the current and next blocks (see Fig. S1/S2 for an example). Each block
$\lfloor\Delta\epsilon\rfloor^{\Gamma}_{I}$ is roughly uniform. The standard
deviation of the mean separation,
$\sigma_{\delta\bar{\omega}}=\sigma_{\delta\omega}/\sqrt{N_{I}}$ ($N_{I}$ is
the number of intervals in the current block), is further calculated to
account for the uncertainty of $\delta\bar{\omega}$. This gives rise to
$T_{I}^{\Gamma}=4\pi/(\delta\bar{\omega}-\sigma_{\delta\bar{\omega}})$. Since
the so-determined $T_{I}^{\Gamma}$ is merely a lower bound, it will be
amplified by a factor of 1.2 for safety, thereby leading to
$T_{I}^{\Gamma}=1.2\times
4\pi/(\delta\bar{\omega}-\sigma_{\delta\bar{\omega}})$.

  3. (c)

Set $T^{\Gamma}=\max_{I}T_{I}^{\Gamma}$.

As a matter of fact, there is no need to go through all the blocks in Step
(c), for the couplings between particle-hole pairs with large energetic
separations are usually very small. Instead, only the first blocks of the
$K$-, $L_{2s}$-, and $L_{2p}$-edge core excitations, partitioned by the
algorithm in Sec. II.4 with $\omega_{G}=0.5\ a.u.$ (see the second column in
Fig. S1/S2), are to be considered.

Finally, given the time series $\\{\mathbf{P}_{\beta}(t_{n})\\}_{n=0}^{N-1}$
($t_{n}=n\Delta t$) of 1PDM, induced by an external field along a direction
$\beta\in\\{x,y,z\\}$, one can calculate many interesting quantities such as
the induced electric dipole moment tensor

$\mu_{\alpha\beta}^{\mathrm{ind}}(t_{n})=-\mathrm{Tr}[\tilde{\mathbf{V}}_{\alpha}\mathbf{P}_{\beta}(t_{n})]+\mathrm{Tr}[\tilde{\mathbf{V}}_{\alpha}\mathbf{P}^{(0)}(0^{-})],\quad
n\in[0,N-1],$ (61)

which can further be transformed to the frequency domain to obtain the
electric dynamical dipole polarizability

$\alpha_{\alpha\beta}(\omega)=\frac{1}{\kappa}\mu_{\alpha\beta}^{\mathrm{ind}}(\omega),$
(62)

and hence the absorption spectrum (dipole strength function)

$S(\omega)=\frac{4\pi\omega}{3c}\Im\mathrm{Tr}\bm{\alpha}(\omega).$ (63)

##  III Pilot Applications

To illustrate the usefulness of the $log_{10}(\Delta
t)$-$log_{10}(|\epsilon_{1s}|)$ relations (see Table 1) for XAS by
CFET4/oCFET4@RT-TDDFT, the HI and PbO molecules are taken here as showcases.
The former is isovalent to the training set (HF, HCl, and HBr) but the latter
is different. The spin-free (sf) partLi  _et al._ (2012, 2014) of the exact
two-component (X2C) HamiltonianKutzelnigg and Liu (2005); Liu and Peng (2009),
the CAM-B3LYPYanai  _et al._ (2004) XC functional, the Dyall-TZP basis
setsDyall (2002) for I and Pb, and the cc-pVTZ basis setsDunning Jr (1989) for
H and O are used in all calculations with the BDF program packageLiu  _et al._
(1997, 2003, 2004); Zhang  _et al._ (2020). The interatomic distance of HI and
PbO are optimized to be 1.602 Å and 1.874 Å, respectively, at the
sf-X2C-CAM-B3LYP level.

To facilitate the dynamics simulation as well as state assignment, the full
molecular symmetry should be used. For HI and PbO, the $A_{1}$ irreducible
representation (irrep) of $C_{\infty v}$ can be used to obtain $A_{1}$ type of
excited states if the applied electric field is along the $z$ (bonding)
direction ($z$-field for short), whereas the $E_{1}$ irrep of $C_{\infty v}$
should be used to obtain $E_{1}$ type of excited states if the applied
electric field is along both the $x$ and $y$ directions. For the latter, it is
operationally much simpler to work with $C_{s}$, the highest subgroup of
$C_{\infty v}$, where $x$ (or $y$) transforms as the totally symmetric irrep
(i.e., $C_{s}(\sigma_{y})$ for $x$ and $C_{s}(\sigma_{x})$ for $y$). Since
$x^{2}$ and $y^{2}$ can be rotated to each other, $\alpha_{xx}(\omega)$ is
strictly identical with $\alpha_{yy}(\omega)$, such that only one of the $x$
and $y$ directions needs to be considered to obtain both $\alpha_{xx}(\omega)$
and $\alpha_{yy}(\omega)$. Nevertheless, only the reflection-symmetric
component of an $E_{1}$ state can be obtained in this case. Note in passing
that it is a must instead of merely a trick to use $C_{s}$ if the applied
$x$\- or $y$-field is very strong. In that case, both $A_{1}$ and $E_{n}$
($n\geq 2$) types of excited states can also be obtained, resulting from the
reflection-symmetric decompositions of $E_{1}^{\beta}\otimes
E_{1}^{\beta}\cdots\otimes E_{1}^{\beta}$, $\beta\in\\{x,y\\}$. In contrast,
the $A_{2}$ type of excited states still cannot be obtained, unless the
angular momentum $l_{z}$ is used in lieu of the electric dipole moment. Having
understood the symmetry issues (i.e., the particle-hole pairs must belong to
the same irreps for both $z$-field under $C_{\infty v}$ and $x/y$-field under
$C_{s}(\sigma_{y/x})$), the time steps $\Delta t$ and total simulation times
$T$ can be determined according to the algorithms described in Secs. II.4 and
II.5, respectively, see Table 3. It is of interest to see that the so-
determined time steps for oCFET4/CFET4 are larger/smaller than the Nyquist
ones ($6.24\times 10^{-5}$ fs for HI and $2.35\times 10^{-5}$ fs for PbO).

Table 3: Predicted time steps ($\Delta t$) and total simulation times ($T$) of CFET4/oCFET4@RT-TDDFT/CAM-B3LYP for core excitations of HI and PbO molecule |  field direction |  $\Delta t$/fs CFET4 |  $\Delta t$/fs oCFET4 |  $T$/fs
---|---|---|---|---
HI | z | $4.38\times 10^{-5}$ | $9.26\times 10^{-5}$ | $4.06$
x | $4.38\times 10^{-5}$ | $9.26\times 10^{-5}$ | $5.28$
PbO | z | $1.66\times 10^{-5}$ | $3.51\times 10^{-5}$ | $6.79$
x | $1.66\times 10^{-5}$ | $3.51\times 10^{-5}$ | $7.05$

The $K$\- and $L$-edge absorption spectra of PbO (HI) can be obtained by the
FD analysisMandelshtam and Taylor (1997) of the CFET4@RT-TDDFT/CAM-B3LYP and
oCFET4@RT-TDDFT/CAM-B3LYP signals for each 1PDM element $P_{ai}(t)$ (more
precisely $P_{ai}(t)+P_{ia}(t)=2\Re P_{ai}(t)$). As can be seen from Figs. 5-7
(Figs. S3-S5), the so-obtained spectra are in excellent agreement with those
by LR-TDDFT/CAM-B3LYP, where the fully symmetrizedPeng  _et al._ (2009)
iterative vector interaction (iVI) approachHuang  _et al._ (2017); Huang and
Liu (2019) is employed to obtain directly the core excitations as interior
roots of LR-TDDFT, without knowing in advance the number or characters of the
roots. To make the comparison even sharper, the $K$-, $L_{2s}$-, and
$L_{2p}$-edge core excitation energies of PbO and HI are further given in
Tables S5-S7 and Tables S8-S10, respectively. It can be seen that the
deviations of CFET4/oCFET4@RT-TDDFT from iVI-LR-TDDFTHuang  _et al._ (2017);
Huang and Liu (2019) are all below 0.05 eV, with the mean absolute error being
0.04 eV for both PbO and HI. This is just the tolerance when training the
$log_{10}(\Delta t)$-$log_{10}(|\epsilon_{1s}|)$ linear relations with respect
to the lowest five bright $K$-edge core excited states of HX (X = F, Cl, Br).
Moreover, it can be seen from Fig. 8 that the predicted time steps are indeed
nearly optimal. As such, the robustness of these relations is fully confirmed.

Having obtained the absorption spectra, it remains to see how individual
particle-hole pairs contribute to the spectra. Instead of usingRepisky  _et
al._ (2015) the imaginary parts of the isotropically averaged induced dipole
moment pairs $\bar{\mu}_{ai}^{\mathrm{ind}}(\omega_{I})$

$\displaystyle\bar{\mu}_{ai}^{\mathrm{ind}}(\omega_{I})$
$\displaystyle=\frac{1}{3}\sum_{\beta\in\\{x,y,z\\}}\tilde{V}^{(1)}_{\beta,ia}d_{\beta,ai}(\omega_{I}),$
(64) $\displaystyle d_{\beta,ai}(\omega_{I})$
$\displaystyle=\int_{-\infty}^{+\infty}[P_{\beta,ai}(t)+P_{\beta,ia}(t)]e^{i\omega_{I}t}dt,$
(65)

the imaginary parts of $\\{d_{\beta,ai}(\omega_{I})\\}$ are employed here to
estimate the compositions $T_{ai}(\omega_{I})$ of state $I$:

$\displaystyle T_{ai}(\omega_{I})$
$\displaystyle=\frac{\bar{P}_{ai}^{2}(\omega_{I})}{\sum_{bj}\bar{P}_{bj}^{2}(\omega_{I})},$
(66) $\displaystyle\bar{P}_{ai}^{2}(\omega_{I})$
$\displaystyle=\frac{1}{3}\sum_{\beta\in\\{x,y,z\\}}[\Im
d_{\beta,ai}(\omega_{I})]^{2}.$ (67)

Noticing that $d_{\beta,ai}(\omega_{I})$ (65) is just the (complex) amplitude
(associated with $\omega_{I}$) in the expression (60) for the signal $2\Re
P_{ai}(t)$ of each particle-hole pair, it can readily be obtained from the FD
analysisMandelshtam and Taylor (1997), which usually requires a much shorter
simulation time than the Fourier transform (65). As can be seen from Tables
S5-S10, the so-obtained transition compositions are indeed very close to those
according to the eigenvectors of iVI-LR-TDDFT. Had
$\Im\bar{\mu}_{ai}^{\mathrm{ind}}(\omega_{I})$ (64) been used, the transition
compositions would be scaled down by the electric dipole moment matrix
elements $\tilde{V}^{(1)}_{ai}$, so as to deteriorate the agreement. As a
matter of fact, it can be seen from equation (104) in Ref. 56 that
$P_{ai}(\omega_{I})$ and $[P_{ai}(-\omega_{I})]^{*}$ ($\equiv
P_{ai}^{*}(\omega_{I})$) are just proportional, respectively, to the upper
($X_{ai,I}$) and lower ($Y_{ai,I}$) components of the $I$th eigenvector of LR-
TDDFT, i.e., $P_{ai}(\omega_{I})\sim
iX_{ai,I}(\mathbf{X}_{I}^{\dagger}\tilde{\mathbf{V}}^{(1)})$ and
$P_{ai}^{*}(\omega_{I})\sim
iY_{ai,I}(\mathbf{Y}_{I}^{\dagger}\tilde{\mathbf{V}}^{(1)})$. Therefore, $\Im
d_{\beta,ai}(\omega_{I})\sim
X_{ai,I}(\mathbf{X}_{I}^{\dagger}\tilde{\mathbf{V}}^{(1)})+Y_{ai,I}(\mathbf{Y}_{I}^{\dagger}\tilde{\mathbf{V}}^{(1)})$
is related directly to the coefficient $\\{C_{ai,I}\\}$ of the excited state
$|I\rangle$ in the basis of singly excited determinants (i.e.,
$|I\rangle=\sum_{ai}C_{ai,I}a_{a}^{\dagger}a_{i}|0\rangle$; cf. equation (113)
in Ref. 56).

Some remarks on the stability and efficiency of the algorithms are now in
order. For the former, it can be seen from Fig. 9(a) that the energy
conservation errors (ECE), $\Delta E_{k}(\Delta t)=E(t_{k})-E(0^{+})$, are
indeed very small for both CFET4@RT-TDDFT/CAM-B3LYP and oCFET4@RT-
TDDFT/CAM-B3LYP, although they grow roughly linearly for each time step as a
result of error accumulation. It can also be seen from Fig. 9(b) that the
ECEs, quantified by $\Delta E_{av}(\Delta t)=\frac{1}{N}\sum_{k=1}^{N}|\Delta
E_{k}(\Delta t)|$, are not monotonic with respect to time step: a larger time
step may lead to a smaller ECE, which is also observed in other
propagtorsGomez Pueyo  _et al._ (2018). It has been attempted to modify the
EPEP2 (EPEP3) predictor-corrector by putting $\mathbf{P}(t_{1}^{[4]})$
($\mathbf{P}(t_{1}^{[6]})$) also into the iteration cycle, by backward Euler
or ETRS propagation from $\mathbf{P}(t_{2}^{[4]})$
($\mathbf{P}(t_{2}^{[6]})$). However, the stability is not improved
discernibly. Anyhow, such small ECEs do not affect the spectral accuracy even
when the total simulation times are doubled (see Table S11).

Since the same total simulation time is used, the relative computational cost
of CFET4@RT-TDDFT and oCFET4@RT-TDDFT can be estimated by the ratio between
their $N_{\mathrm{FBS}}/\Delta t$ values (cf. Tables 2 and 3). In the present
simulations, both EPEP2 in CFET4 and EPEP3 in oCFET4 require two cycles to
converge at each time step, such that the value of $N_{\mathrm{FBS}}/\Delta t$
is equivalent to the number of Fock builds per unit time ($N_{\mathrm{FBT}}$).
It can then readily be checked that CFET4@RT-TDDFT is ca. 1.4 times as
expensive as oCFET4@RT-TDDFT for both HI and PbO. That is, the computational
overhead of oCFET4@RT-TDDFT over CFET4@RT-TDDFT per time step is compensated
for by its larger time step and hence fewer number of steps. As can be seen
from Fig. 10, oCFET4@RT-TDDFT remains to be more efficient than CFET4@RT-TDDFT
for a wide range of spectroscopic accuracy. In contrast, MP4@RT-TDDFT is much
more expensive. As a matter of fact, MP4@RT-TDDFT is not competitive even for
systems containing only light atoms. For instance, it is ca. 2.2 times as
expensive as oCFET4@RT-TDDFT to achieve an accuracy of $<0.05$ eV for the
$K$-edge excitations of HF. Admittedly, even the most efficient scheme
considered here (i.e., oCFET4@RT-TDDFT) is more than four orders of magnitude
more expensive than iVI-LR-TDDFTHuang  _et al._ (2017); Huang and Liu (2019)
(see Table S13), solely due to the necessity of using a very small time step
to fully resolve the core excitation spectra. However, the situation may be
changed for wide absorption spectra of large systems with high densities of
states, as demonstrated beforeTussupbayev  _et al._ (2015).

---
(a)
(b)
(c)
Figure 5: $K$-edge absorption spectra of PbO simulated by (a) CFET4@RT-
TDDFT/CAM-B3LYP, (b) oCFET4@RT-TDDFT/CAM-B3LYP, and (c) LR-TDDFT/CAM-B3LYP.
---
(a)
(b)
(c)
Figure 6: $L_{2s}$-edge absorption spectra of PbO simulated by (a) CFET4@RT-
TDDFT/CAM-B3LYP, (b) oCFET4@RT-TDDFT/CAM-B3LYP, and (c) LR-TDDFT/CAM-B3LYP.
---
(a)
(b)
(c)
Figure 7: $L_{2p}$-edge absorption spectra of PbO simulated by (a) CFET4@RT-
TDDFT/CAM-B3LYP, (b) oCFET4@RT-TDDFT/CAM-B3LYP, and (c) LR-TDDFT/CAM-B3LYP.
---
(a)
(b)
Figure 8: Mean absolute error (MAE) of (a) CFET4/oCFET4@RT-TDDFT/CAM-B3LYP and
(b) MP4@RT-TDDFT/CAM-B3LYP as function of time steps for $K$-edge excitations
of PbO under the z-field. The optimal time steps for MAE $=0.05$ eV are marked
as asterisks, while the predicted time steps (cf. Table 3) are marked as
triangles.
---
(a)
(b)
Figure 9: (a) Energy conservation error $\Delta E_{k}(\Delta
t)=E(t_{k})-E(0^{+})$ as function of simulation time and (b) average energy
conservation error $\Delta E_{av}(\Delta t)=\frac{1}{N}\sum_{k=1}^{N}|\Delta
E_{k}(\Delta t)|$ as function of time steps. CFET4/oCFET4@RT-TDDFT/CAM-B3LYP
results for PbO under the z-field.
---
Figure 10: Number of Fock builds per unit time ($N_{\mathrm{FBT}}$) as
function of the mean absolute error (MAE) by MP4/CFET4/oCFET4@RT-
TDDFT/CAM-B3LYP for $K$-edge excitations of PbO under the z-field (cf. Table
S12).

##  IV Conclusion

A very robust AutoPST approach has been proposed to run relativistic RT-TDDFT
simulations of XAS of systems composed of any elements in the Periodic Table.
It is composed of the following ingredients: (1) the just good enough and most
efficient propagator oCFET4Alvermann and Fehske (2011) along with the EPEP3
predictor-corrector, (2) very robust linear relations and block aliasing
conditions for determining automatically the largest possible time steps, (3)
_a priori_ estimates of the total simulation times, and (4) the use of full
molecular symmetry. Although only scalar relativity has been taken into
account in the pilot applications, AutoPST should work well also for fully
relativistic RT-TDDFT, even beyond the linear response regime.

###### Acknowledgements.

This work was supported by the Key-Area Research and Development Program of
Guangdong Province (Grant No. 2020B0101350001), National Natural Science
Foundation of China (Grant Nos. 21833001 and 21973054), and Mountain Tai
Climbing Program of Shandong Province.

## Data Availability Statement

The data that support the findings of this study are available within the
article and its supplementary material.

## Conflicts of interest

There are no conflicts to declare.

## Supporting Information

Orbital energies, core excitation energies, and absorption spectra of HI and
PbO.

##  Appendix A Algorithms for CFET4 and oCFET4

Algorithm 1 Fourth-order commutator-free exponential time-propagator (CFET4)
with two exponentials (equation (43) in Ref.77)

Subroutine CFET4($\mathbf{P}(0)$, $\mathbf{F}(t_{1})$, $\mathbf{F}(t_{2})$)

1:Inputs:

2: $\mathbf{P}(0)$, $\mathbf{F}(t_{1})$, $\mathbf{F}(t_{2})\leftarrow
t_{1}=(\frac{1}{2}-\frac{\sqrt{3}}{6})\Delta t$,
$t_{2}=(\frac{1}{2}+\frac{\sqrt{3}}{6})\Delta t$

3:$\mathbf{U}_{1}(\Delta t,0)\leftarrow
e^{-i\left(\frac{3+2\sqrt{3}}{12}\mathbf{F}(t_{1})+\frac{3-2\sqrt{3}}{12}\mathbf{F}(t_{2})\right)\Delta
t}$

4:$\mathbf{P}(\Delta t)\leftarrow\mathbf{U}_{1}(\Delta
t,0)\mathbf{P}(0)\mathbf{U}_{1}^{\dagger}(\Delta t,0)$

5:$\mathbf{U}_{2}(\Delta t,0)\leftarrow
e^{-i\left(\frac{3-2\sqrt{3}}{12}\mathbf{F}(t_{1})+\frac{3+2\sqrt{3}}{12}\mathbf{F}(t_{2})\right)\Delta
t}$

6:$\mathbf{P}(\Delta t)\leftarrow\mathbf{U}_{2}(\Delta t,0)\mathbf{P}(\Delta
t)\mathbf{U}_{2}^{\dagger}(\Delta t,0)$

7:Outputs:

8: $\mathbf{P}(\Delta t)$

Algorithm 2 Optimized fourth-order commutator-free exponential time-propagator
(oCFET4) with three exponentials [see equation (43) in Ref.75. The
coefficients $x_{ik}$ in Eq. (34) are determined by expressing the
coefficients $\mathbf{A}_{n}$ defined in equation (25) in Ref. 75 in terms of
the matrix moments $\mathbf{A}_{j}^{[2s]}$ (24)]

Subroutine
oCFET4($\mathbf{P}(0),\mathbf{F}(t_{1}),\mathbf{F}(t_{2}),\mathbf{F}(t_{3})$)

1:Inputs:

2: $\mathbf{P}(0)$ $\mathbf{F}(t_{1})\leftarrow
t_{1}=(\frac{1}{2}-\frac{\sqrt{15}}{10})\Delta t$ $\mathbf{F}(t_{2})\leftarrow
t_{2}=\frac{1}{2}\Delta t$ $\mathbf{F}(t_{3})\leftarrow
t_{3}=(\frac{1}{2}+\frac{\sqrt{15}}{10})\Delta t$

3:$a\leftarrow\left(\begin{array}[]{ccc}\frac{37}{240}-\frac{10\sqrt{15}}{261}&-\frac{1}{30}&\frac{37}{240}+\frac{10\sqrt{15}}{261}\\\
-\frac{11}{360}&\frac{23}{45}&-\frac{11}{360}\\\
\frac{37}{240}+\frac{10\sqrt{15}}{261}&-\frac{1}{30}&\frac{37}{240}-\frac{10\sqrt{15}}{261}\end{array}\right)$

4:$\mathbf{U}_{1}(\Delta t,0)\leftarrow e^{-i\Delta
t(a_{31}\mathbf{F}(t_{1})+a_{32}\mathbf{F}(t_{2})+a_{33}\mathbf{F}(t_{3}))}$

5:$\mathbf{P}(\Delta t)\leftarrow\mathbf{U}_{1}(\Delta
t,0)\mathbf{P}(0)\mathbf{U}_{1}^{\dagger}(\Delta t,0)$

6:$\mathbf{U}_{2}(\Delta t,0)\leftarrow e^{-i\Delta
t(a_{21}\mathbf{F}(t_{1})+a_{22}\mathbf{F}(t_{2})+a_{23}\mathbf{F}(t_{3}))}$

7:$\mathbf{P}(\Delta t)\leftarrow\mathbf{U}_{2}(\Delta t,0)\mathbf{P}(\Delta
t)\mathbf{U}_{2}^{\dagger}(\Delta t,0)$

8:$\mathbf{U}_{3}(\Delta t,0)\leftarrow e^{-i\Delta
t(a_{11}\mathbf{F}(t_{1})+a_{12}\mathbf{F}(t_{2})+a_{13}\mathbf{F}(t_{3}))}$

9:$\mathbf{P}(\Delta t)\leftarrow\mathbf{U}_{3}(\Delta t,0)\mathbf{P}(\Delta
t)\mathbf{U}_{3}^{\dagger}(\Delta t,0)$

10:Outputs:

11: $\mathbf{P}(\Delta t)$

Algorithm 3 Exponential prediction of density matrix and exponential
correction of density matrix at two Gauss-Legendre points (EPEP2) for CFET4

1:Inputs:

2: $\mathbf{P}(0^{-})\leftarrow$ ground state DFT density matrix
$t_{0},t_{f},\Delta t$ $\tilde{\mathbf{V}}^{(1)}\leftarrow$ initial
perturbation

3:Initialize:

4: $N\leftarrow(t_{f}-t_{0})/\Delta t$ $\mathbf{U}(0^{+},0^{-})\leftarrow
e^{-i\tilde{\mathbf{V}}^{(1)}}$
$\mathbf{P}(0)\equiv\mathbf{P}(0^{+})\leftarrow\mathbf{U}(0^{+},0^{-})\mathbf{P}(0^{-})\mathbf{U}^{\dagger}(0^{+},0^{-})$

5:for step $k=0$ to $N-1$ do$\triangleright$ Propagation loop

6: $t_{k}\leftarrow k\Delta t$

7: $t_{a}\leftarrow t_{k}+(\frac{1}{2}-\frac{\sqrt{3}}{6})\Delta t$

8: $t_{b}\leftarrow t_{k}+(\frac{1}{2}+\frac{\sqrt{3}}{6})\Delta t$

9:
$\mathbf{P}_{1}(t_{a})\leftarrow\mathrm{ETRS\\_1}(\mathbf{P}(t_{k}),t_{k},t_{a})$
$\triangleright$ Step 1 (Algorithm 4)

10:
$\mathbf{F}_{1}(t_{a})\leftarrow\mathrm{Fock}(\mathbf{P}_{1}(t_{a}))$$\triangleright$
Step 2

11: $\mathbf{U}_{1}(t_{b},t_{a})\leftarrow
e^{-i\mathbf{F}_{1}(t_{a})(t_{b}-t_{a})}$

12:
$\mathbf{P}_{1}(t_{b})\leftarrow\mathbf{U}_{1}(t_{b},t_{a})\mathbf{P}_{1}(t_{a})\mathbf{U}_{1}^{\dagger}(t_{b},t_{a})$$\triangleright$
Step 3

13:
$\mathbf{F}_{1}(t_{b})\leftarrow\mathrm{Fock}(\mathbf{P}_{1}(t_{b}))$$\triangleright$
Step 4

14:
$\mathbf{P}_{2}(t_{k+1})\leftarrow\mathrm{CFET4}(\mathbf{P}(t_{k}),\mathbf{F}_{1}(t_{a}),\mathbf{F}_{1}(t_{b}))$$\triangleright$
Step 5 (Algorithm 1)

15: for iteration $i=2$ to $N_{\mathrm{it}}$ do $\triangleright$ Loop until
$\|\mathbf{P}_{i+1}(t_{k+1})-\mathbf{P}_{i}(t_{k+1})\|_{F}<M\xi$

16:
$\mathbf{F}_{i}(t_{k+1})\leftarrow\mathrm{Fock}(\mathbf{P}_{i}(t_{k+1}))$$\triangleright$
Step 6

17: $\mathbf{U}_{i}(t_{b},t_{k+1})\leftarrow
e^{-i\mathbf{F}_{i}(t_{k+1})(t_{k+1}-t_{b})}$

18:
$\mathbf{P}_{i}(t_{b})\leftarrow\mathbf{U}_{i}(t_{b},t_{k+1})\mathbf{P}_{i}(t_{k+1})\mathbf{U}_{i}^{\dagger}(t_{b},t_{k+1})$$\triangleright$
Step 7

19:
$\mathbf{F}_{i}(t_{b})\leftarrow\mathrm{Fock}(\mathbf{P}_{i}(t_{b}))$$\triangleright$
Step 4

20:
$\mathbf{P}_{i+1}(t_{k+1})\leftarrow\mathrm{CFET4}(\mathbf{P}(t_{k}),\mathbf{F}_{1}(t_{a}),\mathbf{F}_{i}(t_{b}))$$\triangleright$
Step 5 (Algorithm 1)

21: end for

22: $\mathbf{P}(t_{k+1})\leftarrow\mathbf{P}_{N_{\mathrm{it}}+1}(t_{k+1})$

23:end for

Algorithm 4 One-step enforced time reversal symmetry propagator (ETRS, Eq.
(39))

Subroutine ETRS_1($\mathbf{P}(t_{0})$, $t_{0}$, $t_{1}$)

1:Inputs:

2: $\mathbf{P}(t_{0})$, $t_{0}$, $t_{1}$

3:$\Delta t\leftarrow t_{1}-t_{0}$

4:$\mathbf{F}(t_{0})\leftarrow\mathrm{Fock}(\mathbf{P}(t_{0}))$

5:$\mathbf{U}_{0}(t_{0}+\frac{\Delta t}{2},t_{0})\leftarrow
e^{-i\mathbf{F}(t_{0})\frac{\Delta t}{2}}$

6:$\mathbf{P}_{1}(t_{0}+\frac{\Delta
t}{2})\leftarrow\mathbf{U}_{0}(t_{0}+\frac{\Delta
t}{2},t_{0})\mathbf{P}(t_{0})\mathbf{U}_{0}^{\dagger}(t_{0}+\frac{\Delta
t}{2},t_{0})$ $\triangleright$ Step 1

7:$\mathbf{U}_{1}(t_{1},t_{0})\leftarrow e^{-i\mathbf{F}_{1}(t_{0})\Delta t}$

8:$\mathbf{P}_{2}(t_{1})\leftarrow\mathbf{U}_{1}(t_{1},t_{0})\mathbf{P}(t_{0})\mathbf{U}_{1}^{\dagger}(t_{1},t_{0})$
$\triangleright$ Step 2

9:for iteration $i=2$ to $N_{\mathrm{it}}$ do$\triangleright$ Loop until
$\|\mathbf{P}_{i+1}(t_{k+1})-\mathbf{P}_{i}(t_{k+1})\|_{F}<M\xi$

10:
$\mathbf{F}_{i}(t_{1})\leftarrow\mathrm{Fock}(\mathbf{P}_{i}(t_{1}))$$\triangleright$
Step 3

11: $\mathbf{U}_{i}(t_{1},t_{0}+\frac{\Delta t}{2})\leftarrow
e^{-i\mathbf{F}_{i}(t_{1})\frac{\Delta t}{2}}$

12: $\mathbf{P}_{i+1}(t_{1})\leftarrow\mathbf{U}_{i}(t_{1},t_{0}+\frac{\Delta
t}{2})\mathbf{P}(t_{0}+\frac{\Delta
t}{2})\mathbf{U}_{i}^{\dagger}(t_{1},t_{0}+\frac{\Delta
t}{2})$$\triangleright$ Step 4

13:end for

14:Outputs:

15: $\mathbf{P}_{N_{\mathrm{it}}+1}(t_{1})$

Algorithm 5 Exponential prediction of density matrix and exponential
correction of density matrix at three Gauss-Legendre points (EPEP3) for the
oCFET4 propagator

1:Inputs:

2: $\mathbf{P}(0^{-})\leftarrow$ ground state DFT density matrix
$t_{0},t_{f},\Delta t$ $\tilde{\mathbf{V}}^{(1)}\leftarrow$ initial
perturbation

3:Initialize:

4: $N\leftarrow(t_{f}-t_{0})/\Delta t$ $\mathbf{U}(0^{+},0^{-})\leftarrow
e^{-i\tilde{\mathbf{V}}^{(1)}}$
$\mathbf{P}(0)\equiv\mathbf{P}(0^{+})\leftarrow\mathbf{U}(0^{+},0^{-})\mathbf{P}(0^{-})\mathbf{U}^{\dagger}(0^{+},0^{-})$

5:for step $k=0$ to $N-1$ do$\triangleright$ Propagation loop

6: $t_{k}\leftarrow k\Delta t$

7: $t_{a}\leftarrow t_{k}+(\frac{1}{2}-\frac{\sqrt{15}}{10})\Delta t$

8: $t_{b}\leftarrow t_{k}+\frac{1}{2}\Delta t$

9: $t_{c}\leftarrow t_{k}+(\frac{1}{2}+\frac{\sqrt{15}}{10})\Delta t$

10: $t_{d}\leftarrow\frac{t_{b}+t_{c}}{2}$

11:
$\mathbf{P}_{1}(t_{a})\leftarrow\mathrm{ETRS\\_1}(\mathbf{P}(t_{k}),t_{k},t_{a})$$\triangleright$
Step 1 (Algorithm 4)

12:
$\mathbf{F}_{1}(t_{a})\leftarrow\mathrm{Fock}(\mathbf{P}_{1}(t_{a}))$$\triangleright$
Step 2

13: $\mathbf{U}_{1}(t_{b},t_{a})\leftarrow
e^{-i\mathbf{F}_{1}(t_{a})(t_{b}-t_{a})}$

14:
$\mathbf{P}_{1}(t_{b})\leftarrow\mathbf{U}_{1}(t_{b},t_{a})\mathbf{P}_{1}(t_{a})\mathbf{U}_{1}^{\dagger}(t_{b},t_{a})$$\triangleright$
Step 3

15:
$\mathbf{F}_{1}(t_{b})\leftarrow\mathrm{Fock}(\mathbf{P}_{1}(t_{b}))$$\triangleright$
Step 4

16: $\mathbf{U}_{1}(t_{c},t_{a})\leftarrow
e^{-i\mathbf{F}_{1}(t_{b})(t_{c}-t_{a})}$

17:
$\mathbf{P}_{1}(t_{c})\leftarrow\mathbf{U}_{1}(t_{c},t_{a})\mathbf{P}_{1}(t_{a})\mathbf{U}_{1}^{\dagger}(t_{c},t_{a})$$\triangleright$
Step 5

18:
$\mathbf{F}_{1}(t_{c})\leftarrow\mathrm{Fock}(\mathbf{P}_{1}(t_{c}))$$\triangleright$
Step 6

19:
$\mathbf{P}_{2}(t_{k+1})\leftarrow\mathrm{oCFET4}(\mathbf{P}(t_{k}),\mathbf{F}_{1}(t_{a}),\mathbf{F}_{1}(t_{b}),\mathbf{F}_{1}(t_{c}))$$\triangleright$
Step 7 (Algorithm 2)

20: for iteration $i=2$ to $N_{\mathrm{it}}$ do$\triangleright$ Loop until
$\|\mathbf{P}_{i+1}(t_{k+1})-\mathbf{P}_{i}(t_{k+1})\|_{F}<M\xi$

21:
$\mathbf{F}_{i}(t_{k+1})\leftarrow\mathrm{Fock}(\mathbf{P}_{i}(t_{k+1}))$$\triangleright$
Step 8

22: $\mathbf{U}_{i}(t_{c},t_{k+1})\leftarrow
e^{i\mathbf{F}_{i}(t_{k+1})(t_{k+1}-t_{c})}$

23:
$\mathbf{P}_{i}(t_{c})\leftarrow\mathbf{U}_{i}(t_{c},t_{k+1})\mathbf{P}_{i}(t_{k+1})\mathbf{U}_{i}^{\dagger}(t_{c},t_{k+1})$$\triangleright$
Step 9

24:
$\mathbf{F}_{i}(t_{c})\leftarrow\mathrm{Fock}(\mathbf{P}_{i}(t_{c}))$$\triangleright$
Step 10

25: $\mathbf{U}_{i}(t_{d},t_{c})\leftarrow
e^{i\mathbf{F}_{i}(t_{c})\frac{t_{c}-t_{b}}{2}}$

26:
$\mathbf{P}_{i}(t_{d})\leftarrow\mathbf{U}_{i}(t_{d},t_{c})\mathbf{P}_{i}(t_{c})\mathbf{U}_{i}^{\dagger}(t_{d},t_{c})$$\triangleright$
Step 11

27:
$\mathbf{F}_{i}(t_{d})\leftarrow\mathrm{Fock}(\mathbf{P}_{i}(t_{d}))$$\triangleright$
Step 12

28: $\mathbf{U}_{i}(t_{b},t_{c})\leftarrow
e^{i\mathbf{F}_{i}(t_{d})(t_{c}-t_{b})}$

29:
$\mathbf{P}_{i}(t_{b})\leftarrow\mathbf{U}_{i}(t_{b},t_{c})\mathbf{P}_{i}(t_{c})\mathbf{U}_{i}^{\dagger}(t_{b},t-c)$$\triangleright$
Step 13

30:
$\mathbf{F}_{i}(t_{b})\leftarrow\mathrm{Fock}(\mathbf{P}_{i}(t_{b}))$$\triangleright$
Step 14

31:
$\mathbf{P}_{i+1}(t_{k+1})\leftarrow\mathrm{oCFET4}(\mathbf{P}(t_{k}),\mathbf{F}_{1}(t_{a}),\mathbf{F}_{i}(t_{b}),\mathbf{F}_{i}(t_{c}))$$\triangleright$
Step 7 (Algorithm 2)

32: end for

33: $\mathbf{P}(t_{k+1})\leftarrow\mathbf{P}_{N_{\mathrm{it}}+1}(t_{k+1})$

34:end for

## References

  * Norman and Dreuw (2018) P. Norman and A. Dreuw, Chem. Rev. 118, 7208 (2018).
  * Gao  _et al._ (2016) J. Gao, A. Grofe, H. Ren, and P. Bao, J. Phys. Chem. Lett. 7, 5143 (2016).
  * Grofe  _et al._ (2017) A. Grofe, X. Chen, W. Liu, and J. Gao, J. Phys. Chem. Lett. 8, 4838 (2017).
  * Yang  _et al._ (2019) L. Yang, A. Grofe, J. Reimers, and J. Gao, Chem. Phys. Lett. 736, 136803 (2019).
  * Zhao  _et al._ (2021) R. Zhao, A. Grofe, Z. Wang, P. Bao, X. Chen, W. Liu, and J. Gao, J. Phys. Chem. Lett. 12, 7409 (2021).
  * Gilbert  _et al._ (2008) A. T. Gilbert, N. A. Besley, and P. M. Gill, J. Phys. Chem. A 112, 13164 (2008).
  * Liu  _et al._ (2014) J. Liu, Y. Zhang, and W. Liu, J. Chem. Theory Comput. 10, 2436 (2014).
  * Derricotte and Evangelista (2015) W. D. Derricotte and F. A. Evangelista, Phys. Chem. Phys. Chem. 17, 14360 (2015).
  * Liang  _et al._ (2017) Y. Liang, J. Vinson, S. Pemmaraju, W. S. Drisdell, E. L. Shirley, and D. Prendergast, Phys. Rev. Lett. 118, 096402 (2017).
  * Derricotte and Evangelista (2017) W. D. Derricotte and F. A. Evangelista, J. Chem. Theory Comput. 13, 5984 (2017).
  * Oakley and Klobukowski (2018) M. S. Oakley and M. Klobukowski, J. Electron Spectrosc. Relat. Phenom. 227, 44 (2018).
  * Ehlert and Klamroth (2020) C. Ehlert and T. Klamroth, J. Comput. Chem. 41, 1781 (2020).
  * Carter-Fenk and Herbert (2020) K. Carter-Fenk and J. M. Herbert, J. Chem. Theory Comput. 16, 5067 (2020).
  * Shea  _et al._ (2020) J. A. Shea, E. Gwin, and E. Neuscamman, J. Chem. Theory Comput. 16, 1526 (2020).
  * Hait and Head-Gordon (2020) D. Hait and M. Head-Gordon, J. Chem. Theory Comput. 16, 1699 (2020).
  * Grofe  _et al._ (2020) A. Grofe, R. Zhao, A. Wildman, T. F. Stetina, X. Li, P. Bao, and J. Gao, J. Chem. Theory Comput. 17, 277 (2020).
  * Fronzoni  _et al._ (2004) G. Fronzoni, M. Stener, A. Reduce, and P. Decleva, J. Phys. Chem. A 108, 8467 (2004).
  * Stener  _et al._ (2003) M. Stener, G. Fronzoni, and M. d. de Simone, Chem. Phys. Lett. 373, 115 (2003).
  * Ray  _et al._ (2007) K. Ray, S. DeBeer George, E. I. Solomon, K. Wieghardt, and F. Neese, Chem.–Eur. J. 13, 2783 (2007).
  * Besley and Asmuruf (2010) N. A. Besley and F. A. Asmuruf, Phys. Chem. Chem. Phys. 12, 12024 (2010).
  * Liang  _et al._ (2011) W. Liang, S. A. Fischer, M. J. Frisch, and X. Li, J. Chem. Theory Comput. 7, 3540 (2011).
  * Van Kuiken  _et al._ (2012) B. E. Van Kuiken, N. Huse, H. Cho, M. L. Strader, M. S. Lynch, R. W. Schoenlein, and M. Khalil, J. Phys. Chem. Lett. 3, 1695 (2012).
  * Li and Liu (2010) Z. Li and W. Liu, J. Chem. Phys. 133, 064106 (2010).
  * Li  _et al._ (2011) Z. Li, W. Liu, Y. Zhang, and B. Suo, J. Chem. Phys. 134, 134101 (2011).
  * Li and Liu (2011) Z. Li and W. Liu, J. Chem. Phys. 135, 194106 (2011), (E) 2013, 138, 029904.
  * Van Kuiken  _et al._ (2013) B. E. Van Kuiken, M. Valiev, S. L. Daifuku, C. Bannan, M. L. Strader, H. Cho, N. Huse, R. W. Schoenlein, N. Govind, and M. Khalil, J. Phys. Chem. A 117, 4444 (2013).
  * Verma  _et al._ (2016) P. Verma, W. D. Derricotte, and F. A. Evangelista, J. Chem. Theory Comput. 12, 144 (2016).
  * Peng  _et al._ (2019) R. Peng, A. V. Copan, and A. Y. Sokolov, TJ. Phys. Chem. A 123, 1840 (2019).
  * Nenov  _et al._ (2019) A. Nenov, F. Segatta, A. Bruner, S. Mukamel, and M. Garavelli, J. Chem. Phys. 151, 114110 (2019).
  * Lopata  _et al._ (2012) K. Lopata, B. E. Van Kuiken, M. Khalil, and N. Govind, J. Chem. Theory Comput. 8, 3284 (2012).
  * Lopata and Govind (2013) K. Lopata and N. Govind, J. Chem. Theory Comput. 9, 4939 (2013).
  * Fernando  _et al._ (2015) R. G. Fernando, M. C. Balhoff, and K. Lopata, J. Chem. Theory Comput. 11, 646 (2015).
  * Sissay  _et al._ (2016) A. Sissay, P. Abanador, F. Mauger, M. Gaarde, K. J. Schafer, and K. Lopata, J. Chem. Phys. 145, 094105 (2016).
  * Bruner  _et al._ (2016) A. Bruner, D. LaMaster, and K. Lopata, J. Chem. Theory Comput. 12, 3741 (2016).
  * Neville and Schuurman (2018) S. P. Neville and M. S. Schuurman, J. Chem. Phys. 149, 154111 (2018).
  * Chen and Lopata (2020) M. Chen and K. Lopata, J. Chem. Theory Comput. 16, 4470 (2020).
  * Li  _et al._ (2020) X. Li, N. Govind, C. Isborn, A. E. DePrince III, and K. Lopata, Chem. Rev. 120, 9951 (2020).
  * Yang  _et al._ (2022) M. Yang, A. Sissay, M. Chen, and K. Lopata, J. Chem. Theory Comput. (2022).
  * Verma and Bartlett (2016) P. Verma and R. J. Bartlett, J. Chem. Phys. 145, 034108 (2016).
  * Zapata  _et al._ (2019) F. Zapata, E. Luppi, and J. Toulouse, J. Chem. Phys. 150, 234104 (2019).
  * Gao  _et al._ (2004) J. Gao, W. Liu, B. Song, and C. Liu, J. Chem. Phys. 121, 6658 (2004).
  * Gao  _et al._ (2005) J. Gao, W. Zou, W. Liu, Y. Xiao, D. Peng, B. Song, and C. Liu, J. Chem. Phys. 123, 054102 (2005).
  * Bast  _et al._ (2009) R. Bast, H. J. A. Jensen, and T. Saue, Int. J. Quantum Chem. 109, 2091 (2009).
  * Komorovsky  _et al._ (2019) S. Komorovsky, P. J. Cherry, and M. Repisky, J. Chem. Phys. 151, 184111 (2019).
  * Konecny  _et al._ (2019) L. Konecny, M. Repisky, K. Ruud, and S. Komorovsky, J. Chem. Phys. 151, 194112 (2019).
  * Wang  _et al._ (2005) F. Wang, T. Ziegler, E. van Lenthe, S. van Gisbergen, and E. J. Baerends, J. Chem. Phys. 122, 204103 (2005).
  * Peng  _et al._ (2005) D. Peng, W. Zou, and W. Liu, J. Chem. Phys. 123, 144101 (2005).
  * Xu  _et al._ (2009a) W. Xu, J. Ma, D. Peng, W. Zou, W. Liu, and V. Staemmler, Chem. Phys. 356, 219 (2009a).
  * Xu  _et al._ (2009b) W. Xu, Y. Zhang, and W. Liu, Sci. China Ser. B-Chem. 52, 1945 (2009b).
  * Kühn and Weigend (2013) M. Kühn and F. Weigend, J. Chem. Theory Comput. 9, 5341 (2013).
  * Kühn and Weigend (2014) M. Kühn and F. Weigend, J. Chem. Phys. 141, 224302 (2014).
  * Mori  _et al._ (2014) K. Mori, T. Goumans, E. Van Lenthe, and F. Wang, Phys. Chem. Chem. Phys. 16, 14523 (2014).
  * Egidi  _et al._ (2016) F. Egidi, J. J. Goings, M. J. Frisch, and X. Li, J. Chem. Theory Comput. 12, 3711 (2016).
  * Egidi  _et al._ (2017) F. Egidi, S. Sun, J. J. Goings, G. Scalmani, M. J. Frisch, and X. Li, J. Chem. Theory Comput. 13, 2591 (2017).
  * Stetina  _et al._ (2019) T. F. Stetina, J. M. Kasper, and X. Li, J. Chem. Phys. 150, 234103 (2019).
  * Liu and Xiao (2018) W. Liu and Y. Xiao, Chem. Soc. Rev. 47, 4481 (2018).
  * Repisky  _et al._ (2015) M. Repisky, L. Konecny, M. Kadek, S. Komorovsky, O. L. Malkin, V. G. Malkin, and K. Ruud, J. Chem. Theory Comput. 11, 980 (2015).
  * Konecny  _et al._ (2018) L. Konecny, M. Kadek, S. Komorovsky, K. Ruud, and M. Repisky, J. Chem. Phys. 149, 204104 (2018).
  * De Santis  _et al._ (2020) M. De Santis, L. Storchi, L. Belpassi, H. M. Quiney, and F. Tarantelli, J. Chem. Theory Comput. 16, 2410 (2020).
  * Goings  _et al._ (2016) J. J. Goings, J. M. Kasper, F. Egidi, S. Sun, and X. Li, J. Chem. Phys. 145, 104107 (2016).
  * Konecny  _et al._ (2016) L. Konecny, M. Kadek, S. Komorovsky, O. L. Malkina, K. Ruud, and M. Repisky, J. Chem. Theory Comput. 12, 5823 (2016).
  * Kasper  _et al._ (2018) J. M. Kasper, P. J. Lestrange, T. F. Stetina, and X. Li, J. Chem. Theory Comput. 14, 1998 (2018).
  * Kasper  _et al._ (2020) J. M. Kasper, T. F. Stetina, A. J. Jenkins, and X. Li, Chem. Phys. Rev. 1, 011304 (2020).
  * Darapaneni  _et al._ (2020) P. Darapaneni, A. M. Meyer, M. Sereda, A. Bruner, J. A. Dorman, and K. Lopata, J. Chem. Phys. 153, 054110 (2020).
  * Casida (1995) M. E. Casida, “Time-dependent density functional response theory for molecules,” in  _Recent Advances in Density Functional Methods: (Part I)_ , edited by D. P. Chong (WORLD SCIENTIFIC, 1995) pp. 155–192.
  * Peng  _et al._ (2009) D. Peng, J. Ma, and W. Liu, Int. J. Quantum Chem. 109, 2149 (2009).
  * Huang  _et al._ (2017) C. Huang, W. Liu, Y. Xiao, and M. R. Hoffmann, J. Comput. Chem. 38, 2481 (2017).
  * Huang and Liu (2019) C. Huang and W. Liu, J. Comput. Chem. 40, 1023 (2019).
  * Neville and Schuurman (2019) S. P. Neville and M. S. Schuurman, J. Chem. Phys. 150, 184115 (2019).
  * Herbst and Fransson (2020) M. F. Herbst and T. Fransson, J. Chem. Phys. 153, 054114 (2020).
  * Kosloff (1988) R. Kosloff, J. Phys. Chem. 92, 2087 (1988).
  * Castro  _et al._ (2004) A. Castro, M. A. Marques, and A. Rubio, J. Chem. Phys. 121, 3425 (2004).
  * Gomez Pueyo  _et al._ (2018) A. Gomez Pueyo, M. A. Marques, A. Rubio, and A. Castro, J. Chem. Theory Comput. 14, 3040 (2018).
  * Magnus (1954) W. Magnus, Commun. Pur. Appl. Math. 7, 649 (1954).
  * Alvermann and Fehske (2011) A. Alvermann and H. Fehske, J. Comput. Phys. 230, 5930 (2011).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J.-A. Oteo, and J. Ros, Phys. Rep. 470, 151 (2009).
  * Blanes and Moan (2006) S. Blanes and P. C. Moan, App. Num. Math. 56, 1519 (2006).
  * Liu (2020a) W. Liu, J. Chem. Phys. 152, 180901 (2020a).
  * Liu (2020b) W. Liu, Sci. Sin. Chim. 50, 1672 (2020b).
  * Baer and Neuhauser (2004) R. Baer and D. Neuhauser, J. Chem. Phys. 121, 9803 (2004).
  * Wang  _et al._ (2007) F. Wang, C. Y. Yam, G. Chen, and K. Fan, J. Chem. Phys. 126, 134104 (2007).
  * Cheng  _et al._ (2006) C.-L. Cheng, J. S. Evans, and T. Van Voorhis, Phys. Rev. B 74, 155112 (2006).
  * Zhu and Herbert (2018) Y. Zhu and J. M. Herbert, J. Chem. Phys. 148, 044117 (2018).
  * Nyquist (1928) H. Nyquist, Transactions of the American Institute of Electrical Engineers 47, 617 (1928).
  * Shannon (1948) C. E. Shannon, Bell System Technical Journal 27, 379 (1948).
  * Weigend and Ahlrichs (2005) F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005).
  * Rappoport and Furche (2010) D. Rappoport and F. Furche, J. Chem. Phys. 133, 134105 (2010).
  * Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  * Lee  _et al._ (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  * Becke (1993) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
  * Stephens  _et al._ (1994) P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Phys. Chem. 98, 11623 (1994).
  * Iikura  _et al._ (2001) H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, J. Chem. Phys. 115, 3540 (2001).
  * Yanai  _et al._ (2004) T. Yanai, D. P. Tew, and N. C. Handy, Chem. Phys. Lett. 393, 51 (2004).
  * Wall and Neuhauser (1995) M. R. Wall and D. Neuhauser, J. Chem. Phys. 102, 8011 (1995).
  * Mandelshtam and Taylor (1997) V. A. Mandelshtam and H. S. Taylor, J. Chem. Phys. 107, 6756 (1997).
  * Pang  _et al._ (1998) J. W. Pang, T. Dieckmann, J. Feigon, and D. Neuhauser, J. Chem. Phys. 108, 8360 (1998).
  * Schelter and Kümmel (2018) I. Schelter and S. Kümmel, J. Chem. Theory Comput. 14, 1910 (2018).
  * Li  _et al._ (2012) Z. Li, Y. Xiao, and W. Liu, J. Chem. Phys. 137, 154114 (2012).
  * Li  _et al._ (2014) Z. Li, Y. Xiao, and W. Liu, J. Chem. Phys. 141, 054111 (2014).
  * Kutzelnigg and Liu (2005) W. Kutzelnigg and W. Liu, J. Chem. Phys. 123, 241102 (2005).
  * Liu and Peng (2009) W. Liu and D. Peng, J. Chem. Phys. 131, 031104 (2009).
  * Dyall (2002) K. G. Dyall, Theor. Chem. Acc. 108, 335 (2002).
  * Dunning Jr (1989) T. H. Dunning Jr, J. Chem. Phys. 90, 1007 (1989).
  * Liu  _et al._ (1997) W. Liu, G. Hong, D. Dai, L. Li, and M. Dolg, Theor. Chem. Acc. 96, 75 (1997).
  * Liu  _et al._ (2003) W. Liu, F. Wang, and L. Li, J. Theor. Comput. Chem. 2, 257 (2003).
  * Liu  _et al._ (2004) W. Liu, F. Wang, and L. Li, “Relativistic density functional theory: The BDF program package,” in  _Recent Advances in Relativistic Molecular Theory_ , edited by K. Hirao and Y. Ishikawa (World Scientific, Singapore, 2004) pp. 257–282.
  * Zhang  _et al._ (2020) Y. Zhang, B. Suo, Z. Wang, N. Zhang, Z. Li, Y. Lei, W. Zou, J. Gao, D. Peng, Z. Pu, Y. Xiao, Q. Sun, F. Wang, Y. Ma, X. Wang, Y. Guo, and W. Liu, J. Chem. Phys. 152, 064113 (2020).
  * Tussupbayev  _et al._ (2015) S. Tussupbayev, N. Govind, K. Lopata, and C. J. Cramer, J. Chem. Theory Comput. 11, 1102 (2015).