An Open-Source Framework for $N$-Electron Dynamics: II. Hybrid Density Functional Theory/Configuration Interaction Methodology

TL;DR

This paper introduces an open-source framework combining TDDFT and CIS methods to analyze and visualize correlated many-electron dynamics efficiently, enabling detailed insights into electron behavior in complex molecules.

Contribution

It extends existing tools by integrating hybrid TDDFT/CIS methodology into an open-source Python package for scalable, detailed electron dynamics analysis.

Findings

Good quantitative agreement with higher-level wave function methods.

Demonstrated scalability on medium-sized organic molecules.

Revealed mechanistic details of charge migration via flux densities.

Abstract

In this contribution, we extend our framework for analyzing and visualizing correlated many-electron dynamics to non-variational, highly scalable electronic structure method. Specifically, an explicitly time-dependent electronic wave packet is written as a linear combination of $N$-electron wave functions at the configuration interaction singles (CIS) level, which are obtained from a reference time-dependent density functional theory (TDDFT) calculation. The procedure is implemented in the open-source Python program detCI@ORBKIT, which extends the capabilities of our recently published post-processing toolbox [J. Comput. Chem. 37 (2016) 1511]. From the output of standard quantum chemistry packages using atom-centered Gaussian-type basis functions, the framework exploits the multi-determinental structure of the hybrid TDDFT/CIS wave packet to compute fundamental one-electron quantities such as difference electronic densities, transient electronic flux densities, and transition dipole moments. The hybrid scheme is benchmarked against wave function data for the laser-driven state selective excitation in LiH. It is shown that all features of the electron dynamics are in good quantitative agreement with the higher-level method provided a judicious choice of functional is made. Broadband excitation of a medium-sized organic chromophore further demonstrates the scalability of the method. In addition, the time-dependent flux densities unravel the mechanistic details of the simulated charge migration process at a glance.

Equations

$$i\frac{\partial }{\partial t}|{\Psi }_{\mathrm{el}}\left(t\right)⟩=\left(\hat{H}-\hat{\mu }\cdot \vec{F}(t)\right)|{\Psi }_{\mathrm{el}}\left(t\right)⟩.$$

$$\hat{H}=-\frac{1}{2}\sum _{i=1}^N{\nabla }_i^2+\sum _{i=1}^N\sum _{j>i}^N\frac{1}{r_{ij}}-\sum _{i=1}^N\sum _{A=1}^{N_A}\frac{Z_A}{r_{Ai}}\text{,}$$

$$|{\Psi }_{\mathrm{el}}\left(t\right)⟩=\sum _{\lambda }{B}_{\lambda }\left(t\right)|{\Phi }_{\lambda }⟩.$$

$$|{\Phi }_{\lambda }^{\mathrm{CI}}⟩=C_{\text{ref}}^{\left(\lambda \right)}|{\varphi }_{\text{ref}}⟩+\sum _{ar}{C}_a^{r\left(\lambda \right)}|{\varphi }_a^r⟩+\sum _{abrs}{C}_{ab}^{rs\left(\lambda \right)}|{\varphi }_{ab}^{rs}⟩+\dots .$$

$$i\frac{\partial {\phi }_a(\mathbf{r},t)}{dt}=\left(-\frac{{\nabla }^2}{2}+v_{\text{KS}}(\mathbf{r},t)\right){\phi }_a(\mathbf{r},t).$$

$$v_{\text{KS}}(\mathbf{r},\mathbf{t})=v_{\text{Hartree}}(\mathbf{r},t)+v_{\text{ext}}(\mathbf{r},t)+v_{\text{xc}}(\mathbf{r},t).$$

$$\displaystyle\left[\left(\begin{array}[]{rr}\mathbf{A}&\mathbf{B}\\ \mathbf{B}^{\dagger}&\mathbf{A}^{\dagger}\end{array}\right)-\omega\left(\begin{array}[]{rr}-\mathbf{I}&0\\ 0&\mathbf{I}\end{array}\right)\right]\left(\begin{array}[]{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)=-\left(\begin{array}[]{c}\delta\mathbf{v}\\ \delta\mathbf{v}^{\dagger}\end{array}\right),$$

$$\frac{\partial }{\partial t}\rho \left(\mathbf{r},t\right)=-\vec{\nabla }\cdot \mathbf{j}\left(\mathbf{r},t\right).$$

$$\mathbf{y}\left(\mathbf{r},t\right)$$

$$\hat{\rho }\left(\mathbf{r}\right)$$

$$\hat{j}\left(\mathbf{r}\right)$$

$$⟨\hat{F}⟩(t)$$

$$=$$

$$\displaystyle\Braket{\Phi_{\lambda}\left missing}{\hat{F}}{\Phi_{\nu}}$$

$$\displaystyle\Braket{\varphi_{a}\left missing}{\hat{\rho}}{\varphi_{b}}=\rho_{ab}(\mathbf{r})=\varphi_{a}(\mathbf{r})\varphi_{b}(\mathbf{r}).$$

$$\mathbf{j}\left(\mathbf{r},t\right)$$

$${\mathbf{J}}_{\lambda \nu }\left(\mathbf{r},t\right)$$

$$\displaystyle\mathbf{j}_{ab}=-\frac{i}{2}\Big{(}\varphi_{a}(\mathbf{r})\vec{\nabla}\varphi_{b}(\mathbf{r})-\varphi_{b}(\mathbf{r})\vec{\nabla}\varphi_{a}(\mathbf{r})\Big{)}$$

$${\phi }_a\left(\mathbf{r}\right)=\sum _{A=1}^{N_A}\sum _{i_A=1}^{n_{\mathrm{AO}}(A)}{D}_{i_A}^{(a)}{\chi }_{i_A}\left(\mathbf{r}-{\mathbf{R}}_A\right),$$

$$\displaystyle\Braket{\varphi_{a}\left missing}{\hat{F}}{\varphi_{b}}=\sum_{A,B}^{N_{{A}}}\sum_{i_{A}=1}^{n_{{\rm AO}}{(A)}}\sum_{j_{B}=1}^{n_{{\rm AO}}{(B)}}D^{(a)}_{i_{A}}D^{(b)}_{j_{B}}\Braket{\chi_{i_{A}}\left missing}{\hat{F}}{\chi_{j_{B}}}.$$

$$|{\Psi }_{\mathrm{el}}\left(t\right)⟩=\frac{1}{\sqrt{2}}\left(|{\Psi }_{\mathrm{g}}⟩e^{-i{E}_{\mathrm{g}}t/\hslash }+|{\Psi }_{\mathrm{e}}⟩e^{-i(E_{\mathrm{e}}t/\hslash +\eta )}\right)$$

$$\displaystyle\begin{aligned} \rho\left(\mathbf{r},t\right)=&\frac{1}{2}\Big{(}\rho_{\rm g}\left(\mathbf{r}\right)+\rho_{\rm e}\left(\mathbf{r}\right)\Big{)}+\rho_{\rm ge}\left(\mathbf{r}\right)\cos(\Delta E\,t/\hbar+\eta),\end{aligned}$$

$$\begin{array}{r}\mathbf{j}\left(\mathbf{r},t\right)=\text{mbox}Im\left[{\mathbf{J}}_{\mathrm{ge}}\left(\mathbf{r}\right)\right]\sin (\Delta Et/\hslash +\eta ),\end{array}$$

$$y_{\mathrm{A}}(t)$$

Full text

An Open-Source Framework for $N$-Electron Dynamics: II. Hybrid Density Functional Theory/Configuration Interaction Methodology

Gunter Hermann
, Vincent Pohl22footnotemark: 2 33footnotemark: 3, and Jean Christophe Tremblay22footnotemark: 2 corresponding authorInstitut für Chemie und Biochemie, Freie Universität Berlin, Takustraße 3, 14195 Berlin, GermanyThese authors contributed equally to this work.

1 Abstract

In this contribution, we extend our framework for analyzing and visualizing correlated many-electron dynamics to non-variational, highly scalable electronic structure method. Specifically, an explicitly time-dependent electronic wave packet is written as a linear combination of $N$-electron wave functions at the configuration interaction singles (CIS) level, which are obtained from a reference time-dependent density functional theory (TDDFT) calculation. The procedure is implemented in the open-source Python program detCI@ORBKIT, which extends the capabilities of our recently published post-processing toolbox [J. Comput. Chem. 37 (2016) 1511]. From the output of standard quantum chemistry packages using atom-centered Gaussian-type basis functions, the framework exploits the multi-determinental structure of the hybrid TDDFT/CIS wave packet to compute fundamental one-electron quantities such as difference electronic densities, transient electronic flux densities, and transition dipole moments. The hybrid scheme is benchmarked against wave function data for the laser-driven state selective excitation in LiH. It is shown that all features of the electron dynamics are in good quantitative agreement with the higher-level method provided a judicious choice of functional is made. Broadband excitation of a medium-sized organic chromophore further demonstrates the scalability of the method. In addition, the time-dependent flux densities unravel the mechanistic details of the simulated charge migration process at a glance.

2 Introduction

Unraveling the flow of electrons inside a molecule out of equilibrium is key to understand its reactivity. Since the pioneering laser experiments by Zewail and co-workers[1, 2], the development of new light sources has now granted access to the indirect observation of electron dynamics on its natural timescale. To shed light on the mechanistic details of this attosecond dynamics, accurate theoretical methods are required that capture the subtle details of the transient electronic structure evolution. Various approaches based on explicitly time-dependent density functional theory (TDDFT) and wave function ansatz have been developed over the years and enjoyed mixed degrees of success. While TDDFT appears as more intuitive and scalable, it was shown to suffer from problems for ultrafast dynamics in strong laser fields. On the other hand, the advantages of wave function-based methods in terms of convergence become rapidly compensated by their unfavorable computational cost. Further, the intuitive picture of electrons flowing on a molecular skeleton can become blurred by correlation effects between the $N$ particles.

This contribution is motivated by the need for a robust, scalable wave function method to investigate ultrafast $N$-electron dynamics in systems of large dimension. The method we advocate is based on a combination of linear-response time-dependent functional theory (LR-TDDFT) and configuration interaction singles (CIS) methodologies, as was introduced recently [3, 4, 5]. In principle, the method is similar to the well-established CIS ansatz, with the exception that the energies and the pseudo-CIS eigenvectors are obtained from a reference LR-TDDFT calculation. This allows to improve the energetic properties of the states while keeping a simple electron-particle picture to describe the transient $N$-electron wave packet. This TDDFT/CIS hybrid formalism inherits the qualities of both underlying methods and ensures the $N$-representability of all reduced density matrices, at all times and under all laser conditions. The ensuing $N$-electron dynamics remains marred by the non-intuitive interpretation of quantities beyond the density itself.

Recently, we demonstrated that correlated electron dynamics can be accurately described by means of the electronic flux density operator and derived one-electron properties. We introduced an open-source framework [6] to post-process multi-determinantal configuration interaction wave functions directly from the output of standard quantum chemistry packages. It thus becomes possible to reconstruct the transient $N$-representative one-electron density and current density (flux density) using a library of transition moments calculated from the multi-determinantal configuration interaction wave functions, yielding an intuitive tool for visualizing and analyzing the correlated electron dynamics. A wide variety of established wave function-based methods are covered, ranging from configuration interaction singles to Full CI via restricted active space CI and multi-configuration self-consistent-field methods. It is the purpose of this work to extend the formalism to the TDDFT/CIS hybrid formalism mentioned above, which should retain the qualities of the wave function ansatz and the scalability of DFT-based schemes with respect to the system size.

In the next section, the hybrid TDDFT/CIS methodology is first introduced, followed by the description of the analysis toolset based on the flux density. The application section reports on benchmark calculations on the LiH molecule and the demonstration of the scalability of the scheme by investigation of the broadband excitation in an organic chromophore. The findings are summarized in the conclusion section. Unless otherwise stated, atomic units are used throughout the manuscript ($\hbar=m_{e}=e=4\pi\varepsilon_{0}=1$).

3 Theory

3.1 Hybrid TDDFT/CI Methodology

The evolution of the electronic state of a molecular system obeys the time-dependent Schrödinger equation[7], which can be written in the clamped nuclei approximation

[TABLE]

The interaction of the molecular dipole $\hat{\mu}$ with an external laser field $F(t)$ is treated here semi-classically. For a system consisting of $N$ electrons and $N_{A}$ nuclei, the field-free electronic Hamiltonian reads

[TABLE]

where $\frac{1}{r_{ij}}=\frac{1}{|\vec{r}_{i}-\vec{r}_{j}|}$ is inter-electronic Coulomb repulsion, and $r_{Ai}$ is the distance between the $i$th electron and nucleus $A$ of charge $Z_{A}$. In this work, an electronic wave packet $\Ket{\Psi_{{\rm el}}\left(t\right)}$ satisfying Eq. (1) is expressed as a linear superposition of stationary electronic states $\Ket{\Phi_{\lambda}}$

[TABLE]

Here, $B_{\lambda}\left(t\right)$ are the expansion coefficients of state $\lambda$, which describe the time-evolution of the wave packet. For molecules in strong laser fields, a large number of stationary electronic states is required to offer a proper description of the $N$-electron dynamics. The equations of motion for the coefficients in Eq. (1), associated with the basis set expansion Eq. (3), can be integrated numerically.

In the time-dependent configuration interaction methodology, the stationary electronic states are chosen as linear combinations of excited configuration state functions

[TABLE]

The expansion parameters $C^{\left(\lambda\right)}$ are associated with the formal excitation of a reference configuration, $\Ket{\phi_{\textrm{ref}}}$, from occupied orbitals $\left\{a,b,c\right\}$ to virtual orbitals $\left\{r,s,t\right\}$. Including all possible excitations leads to the exact Full CI limit. The reference and excited configurations are defined as Slater determinants, which builds antisymmetrized products of one-electron spin orbitals $\Ket{\varphi_{a}}$. Note that, in the time-dependent configuration interaction (TDCI) methodology in the form presented above, the field-free electronic Hamiltonian is considered to be diagonal in the basis of CI eigenstates at a given level of theory. The matrix elements of the dipole operator can be computed from the knowledge of these eigenfunctions, which serve as a basis for the variational representation of the molecule-field interaction.

For large molecules, it is customary to truncate the CI expansion to a chosen maximum rank of excitations (e.g., CI Singles or CI Singles Doubles) in order to reduce the number of possible excited configurations. Unfortunately, this often compromises the energetic description of the excited states. To circumvent this limitation while keeping the problem computationally tractable, Sonk and Schlegel [3] first recognized that only excitation energies and transition dipole moments are required to perform TDCI simulations. These can be obtained from linear-response time-dependent density functional theory (LR-TDDFT). To generalize this approach, it was proposed to use the solutions of the LR-TDDFT calculation to generate a basis of pseudo-CI eigenstates [5, 4]. All required information for a TDCI simulation is thus available from the output of standard quantum chemistry programs, provided excitations are performed from the ground state.

According to the Runge-Gross theorem, it is possible to recast the $N$-electron Schrödinger equation and calculate all observables from the sole knowledge of the one-electron density. Using the Kohn-Sham ansatz for the density, the $N$-electron time-dependent Schrödinger equation can be mapped onto a one-electron equation for the orbitals

[TABLE]

The time-dependent Kohn-Sham potential $v_{\textrm{KS}}(\mathbf{r,t})$ contains the classical electrostatic interaction ($v_{\textrm{Hartree}}(\mathbf{r},t)$), an external potential ($v_{\textrm{ext}}(\mathbf{r},t)$), and an exchange-correlation contribution ($v_{\textrm{xc}}(\mathbf{r},t)$), i.e.,

[TABLE]

In explicit TDDFT, the Kohn-Sham potential is usually assumed to be local in time. A celebrated success of TDDFT comes from its linear-response formulation, which allows to accurately compute spectral properties of large molecules. For this endeavor, the response kernel of the electron density to an external weak, long wavelength perturbation can be evaluated from the electric susceptibility of the ground state. The search for the poles of the response function can be recast as an eigenvalue problem of the form

[TABLE]

where $\delta\mathbf{v}$ is the response of the system state to the perturbation, $\delta v_{\textrm{ext}}(\mathbf{r},t)$. The elements of matrices $\mathbf{A}$ and $\mathbf{B}$ are obtained from the orbital energies and integrals over the exchange-correlation kernel, see Eq. (6). At the resonance frequencies $\omega$, where the response vanishes ($\delta\mathbf{v}=0$), the solution of the Casida Eq. (15) yields simultaneously the excitations and de-excitations amplitudes, $\mathbf{X}$ and $\mathbf{Y}$. In the present work, we make use of the fact that these are usually given in the output of standard quantum chemistry programs, together with the excitation energies and the oscillator strengths.

From Eq. (15), it is possible to define pseudo-CI Singles eigenvectors in the Tamm-Dancoff approximation, which consists in neglecting the off-diagonal blocks $\mathbf{B}$. This procedure can alter the quality of the energetic properties of the excited states. On the other hand, the dominant characters present in the pseudo-CI eigenvectors are often not strongly affected by this approximation. In the TDDFT/CI procedure, we thus advocate using directly the transition energies and amplitudes obtained from a LR-TDDFT calculation to take advantage of the full solution of Eq. (15) and to obtain a good energetic description of the excited states. A separate Tamm-Dancoff calculation may be used to confirm the character of the excited states. The LR-TDDFT excitation amplitudes are then re-orthonormalized using a modified Gram-Schmidt procedure to define a pseudo-CI basis for the TDCIS dynamics. All properties not directly deriving from the energies can be subsequently calculated at the CIS level of theory using the orbitals and the pseudo-CI eigenstates, which are treated as configuration interaction singles expansions. Note that the Slater determinants $\left\{\phi_{\textrm{ref}},\phi_{a}^{r}\right\}$ are constructed from Kohn-Sham orbitals. Importantly, all the information required to reconstruct these KS-orbitals and the $N$-electron pseudo-eigenfunctions are directly accessible from the output of standard quantum chemistry packages. As a consequence, only the evaluation of one-electron integrals is required to generate a library of molecular properties and transition moments of various one-electron operators, which can be used to characterize the properties of transient wave packets, as explained below.

3.2 Analysis Tools for Electron Dynamics

For the analysis of the $N$-electron dynamics, we propose using a set of tools composed from the one-electron density, $\rho\left(\mathbf{r},t\right)$, and the associated electronic flux density, $\mathbf{j}\left(\mathbf{r},t\right)$. These are related by the electronic continuity equation

[TABLE]

Whereas the electron density gives information about the probability distribution of the electron, the flux density yields complementary information about the phase of the electronic wave packet. This in turn reveals the mechanistic aspects of the time-evolution of the one-electron density. The one-electron density can be used to define the electron flow, $\frac{\partial}{\partial t}\rho\left(\mathbf{r},t\right)$, as the left-hand-side of the continuity equation. The difference density, $\mathbf{y}\left(\mathbf{r},t\right)$, is a widespread quantity used for visualization purposes, and it can be obtained by integrating the electron flow from a chosen initial condition $\rho\left(\mathbf{r},0\right)$, i.e.,

[TABLE]

We will resort to both quantities in later analyses.

In operator form, the one-electron density and the electronic flux density respectively read

[TABLE]

where $\mathbf{r}$ is an observation point, $\delta\left(\mathbf{r}-\mathbf{r}_{k}\right)=\delta_{k}(\mathbf{r})$ is the Dirac delta distribution at the position $\mathbf{r}_{k}$ of electron $k$, and $\hat{p}_{k}=-i\vec{\nabla}_{k}$ is the associated momentum operator. In general, the expectation value of any one-electron operator $\hat{F}$ can be expressed using Eqs. (3) and (4) as

[TABLE]

Evaluation of the matrix elements \Braket{\Phi_{\lambda}\left missing}{\hat{F}}{\Phi_{\nu}} can be done by exploiting the structure of the functions $\left\{\Phi_{\lambda},\Phi_{\nu}\right\}$. In the hybrid TDDFT/CIS methodology, these take the form of singly excited configurations, i.e., the truncation of Eq. (4) at the singles level. The matrix elements in the basis of singly excited configurations read

[TABLE]

where $\overline{a}\in\{1,2,\dots,a-1,r,a+1,\dots\}$ denote the occupied spin orbitals of the configuration state function $\Ket{\phi_{a}^{r}}$. Note that we make use of the Slater-Condon rules[8, 9, 10] to resolve Eq. (23) in terms of one-electron integrals in the basis of the spin orbitals, $\varphi_{a}(\mathbf{r})$. The transition moments between spin orbitals are usually computed in the spin-free representation by first integrating over the spin coordinates. Specifically, the expectation value for the electron density requires the following integrals

[TABLE]

The electronic flux density for a wave packet of the form Eq. (3) can be formulated as

[TABLE]

which can be calculated by exploiting the CIS structure of the eigenfunctions. The transition electronic flux density from state $\lambda$ to state $\nu$ is denoted $\mathbf{J}_{\lambda\nu}\left(\mathbf{r},t\right)$, which simplifies using the Slater-Condon rules to

[TABLE]

where

[TABLE]

are molecular orbital (MO) electronic transition flux densities from MO $\varphi_{a}(\mathbf{r})$ to MO $\varphi_{b}(\mathbf{r})$.

As one of the most widespread bases used in quantum chemistry, we specialize here to spatial MO defined as linear combination of atom-centered orbitals (MO-LCAO for “Molecular Orbital - Linear Combination of Atomic Orbitals”)

[TABLE]

where $D^{(a)}_{i_{A}}$ is the $i_{A}$th expansion coefficient for MO $a$. The atomic orbitals $\chi_{i_{A}}$ are expressed as a function of the Cartesian coordinates of one electron $\mathbf{r}$ and the spatial coordinates $\mathbf{R}_{A}$ of nucleus $A$. $N_{{A}}$ labels the number of atoms and $n_{{\rm AO}}{(A)}$ is the number of atomic orbitals on atom $A$. Using the MO-LCAO ansatz, the transition moments between spin orbitals read

[TABLE]

The MO-LCAO coefficients $D^{(a)}_{iA}$ and the definition of the atomic orbitals can be read directly from the output of standard quantum chemistry program packages. All required derivatives and integrals in the atomic orbital basis are computed analytically using our Python post-processing toolbox ORBKIT[11], with which the molecular orbital density (cf. Eq. (24)) and the molecular orbital electronic flux density (cf. Eq. (27)) can then be projected on an arbitrary grid. Combining the information in this list with the occupation patterns of the quasi-CI eigenvectors associated with the excited states obtained at the LR-TDDFT level of theory, it is possible to create a library of transition moments between CI-states to be used in the dynamics. Note that the transition dipole moments are also computed using the same information and exploiting the multi-determinantal structure of the $N$-electron basis functions, cf. Eqs. (3) and (4). The analysis tools for the hybrid TDDFT/CIS methodology are implemented, along with various other one-electron quantities, in a recently introduced open-source Python framework detCI@ORBKIT, available at https://github.com/orbkit/orbkit/. The program requires a preliminary LR-TDDFT calculation using Gaussian-type atom-centered orbitals. There is no restriction for the choice of functional. Currently, the code supports the GAMESS[12] and TURBOMOLE[13] formats. Our program then computes matrix elements of one-electron operators, projects them on an arbitrary grid, and stores them in a library to be used for analyzing the $N$-electron dynamics. The framework detCI@ORBKIT is written in Python, simplifying its portability on different platforms and offering efficient standard libraries for visualization purposes. Implementation details are given elsewhere [14]. The dynamics program is not part of the standard implementation and can be performed using either a user-written code or, e.g., the Matlab WAVEPACKET package[15]. In the present work, all dynamical simulations were performed using GLOCT, an in-house implementation of a propagator for the reduced-density matrix and related quantities [16, 17].

4 Results and Discussion

To demonstrate the capabilities of detCI@ORBKIT, we perform the analysis of correlated electron dynamics in two selected molecular systems. First, the charge transfer process in lithium hydride is studied to benchmark the quality of the TDDFT/CIS description against Full CI results. The electron migration in an alizarin dye induced by broadband laser excitation is then used as an example to demonstrate the scalability of the method.

4.1 Benchmark: Charge Transfer in LiH

In lithium hydride, charge migration can be initiated, e.g., by laser excitation from the molecular ground state $\Ket{\Psi_{\rm g}}$ to the first excited state $\Ket{\Psi_{\rm e}}$. The charge transfer mechanism can be understood from the analysis of the superposition state

[TABLE]

which leads to the time-dependent one-electron density

[TABLE]

where $\Delta E=E_{{\rm e}}-E_{{\rm g}}$, and $\eta$ is the relative phase. It is set to $\eta=\pi$ in the present example. Similarly, the electronic flux density takes the form

[TABLE]

where the transition electronic flux density $\mathbf{J}_{\rm ge}\left(\mathbf{r}\right)$ is obtained from Eq. (26). The charge transfer dynamics can be thus rationalized in terms of the static transition moment between the two states involved. The electron densities, $\rho_{\rm g}\left(\mathbf{r}\right)$ and $\rho_{\rm e}\left(\mathbf{r}\right)$, of the ground state X${}^{1}\Sigma^{+}$ and the charge transfer state A${}^{1}\Sigma^{+}$ and the transition density between both, $\rho_{\rm ge}\left(\mathbf{r}\right)$, are computed by combining the MO contributions obtained from Eq. (24). The Kohn-Sham orbitals, the pseudo-CI eigenvectors, and the associated LR-TDDFT excitation energies are computed using an aug-cc-pVTZ at the B2-PLYP level of theory, as implemented in TURBOMOLE.[13] The character of the charge transfer state is found to be dominated by the HOMO-LUMO transition (see Fig. 1 (right side)). This is in good agreement with the character determined from CIS and Full CI calculations, both performed with PSI4[18] using the identical basis set. The corresponding frontier orbitals from the Hartree-Fock reference are shown in the left side of Fig. 1. It can be seen that the HOMO is similar in both cases, while the LUMO is more delocalized at the B2-PLYP level of theory. We will show below that this difference has only a marginal influence on the electron dynamics.

A great advantage of LR-TDDFT over CIS is the improved energetic description of the excited states at virtually the same computational cost. This is demonstrated by the good agreement of the excitation energies at the B2-PLYP level of theory ($\Delta E=3.51$ eV) with the Full CI reference ($\Delta E=3.56$ eV), as compared to the rather poor value for the truncated CIS expansion ($\Delta E=4.04$ eV). Since the excitation is dominated by single excitations in all cases, the transition energies will affect mostly the timescale of the dynamics. Further, considering the similarities between the MOs involved in the dominant transition, a similar dynamical behavior is to be expected. This is indeed what is observed in Fig. 2, where the left panels report the flux density $\mathbf{j}\left(\mathbf{r},t\right)$ at time $t=\tau/4$ computed using Eq. (32) and the right panels show the difference density $\mathbf{y}\left(\mathbf{r},t\right)$ computed from the one-electron density, Eq. (31). In the top right panel, the benchmark Full CI calculation reveals that the charge is transferred from the hydrogen atom (electron density depletion region in blue) to the lithium ion (red region). Some $p$-like regions of increasing density can also be recognized around the lithium atom. The same features are also observed for both the CIS (central panels) and hybrid (lower panels) method, while the magnitude of the difference density is larger around the hydrogen for the two single electron approaches.

The electronic flux density maps are depicted as streamlines in the left panels of Fig. 2, where the blue shades indicate its magnitude. The charge is seen to flow from the polarizable hydrogen around and towards the back of the hard lithium. The same large vortex around the lithium ion is observed for the Full CI benchmark, the CIS, and the hybrid TDDFT/CIS approach. The main quantitative differences between the methods are located in the low-density regions, e.g., at $x>5~{}a_{0}$, where the Full CI benchmark predicts a flux almost parallel to the molecular axis. The critical point (at $x\sim 2~{}a_{0}$) between the lithium and hydrogen atoms also appears to be slightly shifted to the right at the CIS and TDDFT/CIS level of theory. In general, it can be said that both single electron excitation ansatz yield a very similar picture of the electron difference density and the associated flux density. This conclusion is likely to hold for all dynamical processes involving $N$-electron eigenstates dominated by a single excitation character.

4.2 Scalability: Electron Migration in Alizarin

In this second example, we demonstrate the computational scalability of the hybrid TDDFT/CIS approach to analyze the correlated electron dynamics for more extended molecular systems. The necessity of such a method is due to the fact that high-level electronic structure theory methods, such as MCSCF, are often not applicable for larger molecules. The CI scheme truncated at the singles excitation represents a simple, intuitive, and computationally cheap approach to compute qualitatively correct excited electronic states.[19, 20] However, it yields inaccurate vertical excitation energies from the ground state.[21] As advocated in the theory section above, a suitable alternative is LR-TDDFT[22], which usually provides better energetic description than CIS while retaining the same quality for the wave function. It can be inferred from the example in the previous subsection, that this will provide an adequate description for a large number of photochemical processes dominated by a single excitation character. In addition, it benefits from the versatility and continuous improvement of density functionals. Proper treatment of the excited states strongly depends on the appropriate choice of a functional, which can be chosen to correctly describe electronic correlations, the dispersive nature, or the charge-transfer character of a given excitation. [23, 24, 25, 26, 27, 28, 29] Fortunately, extensive experience has been accumulated over the years concerning the applicatibility of each functional in specific chemical contexts. For example, it is known that non-local exchange improves greatly the description of charge transfer states [20].

To show the scalability of the hybrid TDDFT/CIS formalism for the analysis of correlated electron dynamics, we initiate a photoinduced ultrafast charge migration process in alizarin. This organic chromophore is used as a $\pi$-conjugated photosensitizer in dye-sensitized solar cells. Prior to the dynamical simulation, the electronic and optical properties of alizarin are determined by means of LR-TDDFT. Therefore, we perform a TURBOMOLE[13] calculation with the B3LYP[30] hybrid functional and the def2-SVP basis set[31, 32] at the equilibrium geometry of alizarin. This setup has been previously proven to yield accurate results for the electronic and spectroscopic properties of such systems.[33, 34, 35, 4] The computed absorption spectrum is depicted in Fig. 3(a) along with the experimentally observed absorption band of free alizarin (dashed black line). The good agreement for the first absorption band between theory (437 nm) and experiment (431 nm) underlines the suitability of TDDFT to model the electronic spectra of medium-sized organic molecules dominated by single excitation character. It is important to recognize that the second absorption band is composed of a multitude of excited states in the UV/VIS range.

In order to simulate an ultrafast charge migration process in alizarin, we proceed to the broadband excitation of all excited states in the energy range between 200 nm and 500 nm (cf. Fig. 3(a) yellow filling). For the promotion of these states from the ground state, a superposition of state-to-state sin2-shaped pulses with a duration 19 fs is constructed. The pulse is adjusted to the parameters of a realistic experimental laser field used in similar investigations to initiate, e.g., ultrafast photoinduced processes in alizarin-TiO2 solar cells.[36, 37] The resulting electric field is shown in the inset of Fig. 3(b). The laser excitation is followed by a 20 fs field-free propagation. The time-evolution of the $N$-electron wave packet (cf. Eqs. (1) and (3)) is accomplished using an adaptive Runge-Kutta algorithm in the interaction picture. The methodology and implementation details are described elsewhere.[38, 16, 39] Fig. 3(b) shows the evolution of the state populations and the applied laser field in the inset. As it is often the case for molecules in strong fields, the population dynamics is very intricate while the laser is on, in part due to important polarization effects and in part due to the number of states that are excited by the broadband laser. To account for the electronic response of the system to the laser field, 25 eigenstates are incorporated in the simulation. After the laser excitation, only twelve states are significantly populated ($P_{\lambda}>0.01$).

To unravel the mechanistic pathways and give an intuitive picture of the electron dynamics, we advocate using the time-dependent electron density, electron flux density, and electronic flow, which are reconstructed from the $N$-electron wave packet in the pseudo-CI eigenvector basis. This can become computationally tedious, since the number of Slater determinants in the wave function expansion increases with the number of occupied and unoccupied orbitals. For alizarin in the current basis set, the $62$ occupied MOs times $248$ virtual MOs correspond to $15\,376$ determinants for each of the excited states. Recalling that reconstruction of the flux density, Eq. (26), requires combining one-electron integrals of all orbitals pairwise, for each pair of eigenstates, this amounts to a tremendous computational task. However, inclusion of the complete set of determinants is not necessary, since the expansion coefficient of many of these is either zero or negligibly small. Moreover, the pseudo-CI eigenvectors are usually dominated by a few determinant contributions. Applying the maximum correspondence principle between determinants and exploiting the Slater-Condon rules lead to significant numerical savings, which are automatized in our implementation and can be controlled by a user-defined convergence threshold.

In order to examine the influence of truncating the complete set of determinants to a physically meaningful subset in the wavefunction ansatz, we choose to use a threshold for the expansion coefficients of the quasi-CI wave functions, Eq. (4). All determinants with a contribution under this value are neglected in the evaluation of the one-electron density and the associated flux density. We define three different thresholds, $\left|C_{a}^{r,(\lambda)}\right|>\{10^{-3},\,10^{-2},\,10^{-1}\}$. To illustrate the computational savings that can be expected, the thirteenth excited state is selected as an example, since it is the most populated state after the broadband laser excitation. The three chosen threshold values retain 926, 31, and 7 Slater determinants in the wave function expansion, respectively. For consistency check, the pruned expansions recover 99.94 $\%$, 99.31 $\%$, and 95.92 $\%$ of the norm of the thirteenth excited state, respectively. These numbers are similar for other excited states. After excluding the contributions lying under the respective threshold, the remaining coefficients are renormalized using a modified Gram-Schmidt algorithm. Regarding the computational effort, the reduction of determinants means a drastic decrease of computational steps. For example, the calculation of the transition electronic flux density $\mathbf{j}_{5\,13}\left(\mathbf{r}\right)$ between the fifth and thirteenth excited state requires $N_{\rm SD}^{5}\times N_{\rm SD}^{13}=573\times 926=530\,598$ determinant combinations for the more strict threshold of $\left|C_{a}^{r,\lambda}\right|>10^{-3}$. This number reduces to $N_{\rm SD}^{5}\times N_{\rm SD}^{13}=3\times 7=21$ determinant combinations for the lenient threshold $\left|C_{a}^{r,\lambda}\right|>10^{-1}$. This simple constraint thus confers great scalability to the method presented here.

To assess the quality of this approximation, the electronic flux densities $\mathbf{j}\left(\mathbf{r},t\right)$ reconstructed using the different thresholds are illustrated in Fig. 4(a),(c),(d) at a characteristic point in time after the laser-pulse excitation (here, $t=25.3$ fs). This analysis could be performed during the laser-field at the zeros of the pulse function to avoid simply representing the contribution of the electric field to the flux density. On the other hand, after the pulse, the flux is solely caused by the coherences between the electronic states, simplifying its analysis. In Fig. 4(b), the corresponding electron flow, $\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$, is displayed additionally for the tight threshold, $\left|C_{a}^{r,(\lambda)}\right|>10^{-3}$, to facilitate the interpretation of the flux densities. For all three wave packet expansions, $\mathbf{j}\left(\mathbf{r},t\right)$ shows nearly identical qualitative features. These include: (1) an electron flow along the bonds mediated by the $\pi$-system, (2) an anti-clockwise flux at the outer right ring, and (3) a charge migration from the outer right ring and the top part of the left ring to the lower hydroxyl group. Due to the cyclic nature of the field-free evolution of the $N$-electron wave packet, the features (2)-(3) exhibits a Rabi-type change of direction during the dynamics. Since the qualitative picture is the same even at all levels of analysis, the very important numerical savings at the crudest level of approximation confer great scalability to the proposed methodology to analyze $N$-electron dynamics in large molecules. Despite this success, one striking difference can be noticed, i.e., the electron migration from the lower hydroxyl group to the neighboring carboxyl group is not fully reproduced with the two larger thresholds (cf. Fig. 4(c) and (d)). This corresponds to a through-space charge transfer and is mediated by a large number of small contributions in the pseudo-CI eigenstate basis. While the major phenomenological characteristics of the correlated electron dynamics are still captured, some minor mechanistic information is lost by reducing the wave functions to their dominant determinant contributions.

To understand the electron dynamics, we extend its analysis to the time-dependent electronic yield. It is defined as the difference between the electron density at a given time and the electron density at $t=0\,{\rm fs}$, integrated over a given volume,

[TABLE]

In Fig. 5, the difference densities projected on right and left aromatic ring of alizarin are reported for times after the laser has been switched off. These reveal intricate synchronous fluctuations of electron density between both rings. The slight asymmetry of the electron redistribution correlates with the asymmetric substituents on the two outer rings: the fluctuation amplitudes in the left ring are larger than in the unsubstituted, rightmost aromatic ring. In Fig. 6, the electronic flux density is plotted for selected characteristic points in order to unveil the mechanism of the charge migration between both rings. The times associated with these snapshots are marked as vertical gray lines in Fig. 5. In panel (a), where the electronic yield is largest on the rightmost aromatic ring ($t=$24.5 fs), the electronic flux density is seen to rotate clockwise, which will create a temporary local magnetic field. In the next snapshot (panel (b), $t=$25.3 fs), the yields in the left and right aromatic rings are equivalent. The electrons are seen to flow laminarly from the right to left aromatic moieties. Starting from the lowest carbon atom of the right ring, electrons move anticlockwise along the bonds of the right ring, passing by the bottom carbonyl group to bottom hydroxyl group below the left aromatic ring, and back to the central carbonyl unit via a through space mechanism. This is a strong indication of the presence of a hydrogen bond, which will lead to a rapid tautomeric hydrogen transfer. At the third time step ($t=$25.9 fs), the electronic yield is maximal in the right ring and the electron flux density rotates clockwise. Contrary to panel (a), the electron flow is not as strongly localized inside the ring, probably due to the asymmetry of the ring substituents. In the last snapshot ($t=$26.4 fs), the density migrates back from the leftmost to the rightmost aromatic ring along a different path: the electrons retain a clockwise orientation in the left ring, mostly flowing from the bottom hydroxyl to the top central carbonyl group. A marginal amount flows from the bottom left hydroxy group over the hydrogen bond to the bottom carbonyl group and into rightmost aromatic ring. For the complete picture of the dynamics after the laser excitation, a short film of the charge migration process is made available online in the Supporting Information.

5 Conclusions

In this paper, we have introduced a novel procedure to analyze and visualize many-electron dynamics from a hybrid time-dependent density functional theory (TDDFT)/ configuration interaction singles (CIS) formalism. The method resorts to a linear-response TDDFT calculation to generate a basis of pseudo-CI eigenvectors and associated energies, which are then used as a basis to describe an $N$-electron dynamics at the CIS level of theory. The time-dependent CIS wave function retains the simple character of coupled electron-hole pairs, which facilitates its interpretation in terms of configuration states while keeping the size of the basis relatively small. This renders the hybrid method amenable to large systems – in fact, any system that can be tackled and described accurately using linear-response TDDFT.

From the evolution of a TDDFT/CIS wave packet, we showed how it is possible to rationalize the $N$-electron dynamics of a system in terms of transition moments of various one-electron operators. These include difference electronic densities, transient electronic flux density maps, and transition dipole moments, which we have implemented in a Python toolbox for postprocessing multi-determinantal wave function data. This new module of our open-source project ORBKIT creates a library of transition moments of user-specified one-electron operators, which are projected on an arbitrary grid. The required information (molecular orbitals, structure, Gaussian atomic basis, pseudo-CI coefficients) can be directly extracted from the output of various quantum chemistry program packages. These are then used to reconstruct quantities that help understanding the flow of electrons in molecules out of equilibrium.

We first applied the hybrid TDDFT/CIS scheme to a test system, the charge transfer in lithium hydride, in order to benchmark the quality of the ansatz against standard CIS and Full CI reference simulations. The results demonstrated that a good choice of the functional improves mostly the energetic description of the charge transfer state, which can be brought close to the Full CI benchmark. On the other hand, the pseudo-CI basis retains a single excitation character, which is found to be similar to the standard CIS reference. Both ansatz agree semi-quantitatively with the higher level wave function description, with the discrepancies mostly found in the regions of low density.

In a second example, the application to the broadband excitation of a prototypical chromophore for dye-sensitized solar cells demonstrated the scalability of the method and the versatility of the new toolkit. In particular, it was found that the main features of the electron flow mechanism can be recovered using a stringent basis pruning strategy, in which each pseudo-CI eigenstate is represented using only a few dominant configurations. By doing so, marginal features of the many-electron dynamics involving, e.g., electron flow through hydrogen bonds, may be lost. Because of the favorable scaling of the method even at tighter convergence Thresholds, it is expected to be applicable to a large number of medium-sized molecules. This will help to understand electron migration processes in various fields, such as nano electronics and light-harvesting applications.

6 Acknowledgment

The authors gratefully thank the Scientific Computing Services Unit of the Zentraleinrichtung für Datenverarbeitung at Freie Universtät Berlin for allocation of computer time, and Hans-Christian Hege for providing the ZIBAmira visualization program. The funding of the Deutsche Forschungsgemeinschaft (project TR1109/2-1) and of the Elsa-Neumann foundation of the Land Berlin are also acknowledged.

7 Keywords

Correlated Electron Dynamics, Time-dependent Density Functional Theory, Electronic Flux Density, Electron Density, Electronic Current Density