An Open-Source Framework for Analyzing $N$-Electron Dynamics: I. Multi-Determinantal Wave Functions

TL;DR

This paper introduces an open-source Python framework for analyzing and visualizing correlated many-electron dynamics in molecules, focusing on electronic flux density derived from multi-determinant wave functions, applicable to charge migration studies.

Contribution

It provides a novel post-processing tool that extracts electronic flux data directly from standard quantum chemistry outputs, enabling detailed analysis of ultrafast electron dynamics.

Findings

Demonstrates broad applicability to charge migration in various molecules.

Assesses convergence of dynamics with respect to theory level and basis set.

Shows robustness of dynamical features across different computational parameters.

Abstract

The aim of the present contribution is to provide a framework for analyzing and visualizing the correlated many-electron dynamics of molecular systems, where an explicitly time-dependent electronic wave packet is represented as a linear combination of $N$-electron wave functions. The central quantity of interest is the electronic flux density, which contains all information about the transient electronic density, the associated phase, and their temporal evolution. It is computed from the associated one-electron operator by reducing the multi-determinantal, many-electron wave packet using the Slater-Condon rules. Here, we introduce a general tool for post-processing multi-determinant configuration-interaction wave functions obtained at various levels of theory. It is tailored to extract directly the data from the output of standard quantum chemistry packages using atom-centered Gaussian-type basis functions. The procedure is implemented in the open-source Python program detCI@ORBKIT, which shares and builds upon the modular design of our recently published post-processing toolbox [J. Comput. Chem. 37 (2016) 1511]. The new procedure is applied to ultrafast charge migration processes in different molecular systems, demonstrating its broad applicability. Convergence of the $N$-electron dynamics with respect to the electronic structure theory level and basis set size is investigated. This provides an assessment of the robustness of qualitative and quantitative statements that can be made concerning dynamical features observed in charge migration simulations.

Figures

Tables

Table 1. Table 1: Excitation energies Δ ​ E Δ 𝐸 \Delta E and transition dipole moments μ r subscript 𝜇 𝑟 \mu_{r} in length gauge for the trihydrogen cation H + 3 superscript subscript absent 3 {}_{3}^{+} for the transition between state 1 ​ A 1 ′ 1 1 superscript superscript subscript 𝐴 1 ′ 1 1{}^{1}A_{1}^{\prime} and state 1 ​ E ′ 1 1 superscript superscript 𝐸 ′ 1 1{}^{1}E^{\prime} at the Full CI of theory for a selection of basis sets. The dipole moment in length gauge μ r subscript 𝜇 𝑟 \mu_{r} is compared to ( μ v ) r subscript subscript 𝜇 𝑣 𝑟 (\mu_{v})_{r} , the dipole moment in length gauge calculated from the dipole moment in velocity gauge (cf. Eq. ( 34 )). Additionally, the number of basis functions is listed for each basis set.

Basis Set	${\mu }_r$ [D]	${\left({\mu }_v\right)}_r$ [D]	$\Delta E$ [eV]	$N_{\mathrm{AO}}$	$N_{\mathrm{DET},g}$	$N_{\mathrm{DET},e}$
STO-3G	2.9336	1.8433	23.6135	003	0003	00000004/4
cc-pVDZ	2.7461	2.7005	19.4211	015	0045	00000068/68
cc-pVTZ	2.7691	2.7545	19.3612	042	0300	0000464/464
cc-pVQZ	2.7668	2.7654	19.3414	090	1292	002011/2000
cc-pV5Z	2.7669	2.7669	19.3354	165	4161	006507/6518
aug-cc-pVDZ	2.7674	2.7265	19.3224	027	0139	0000212/212
aug-cc-pVTZ	2.7683	2.7638	19.3223	069	0817	001248/1222
aug-cc-pVQZ	2.7687	2.7675	19.3256	138	3046	004643/4626
aug-cc-pV5Z	2.7668	2.7675	19.3277	240	9448	13836/13680
Ref. [89]			19.3289	600

Equations

$$i\frac{\partial }{\partial t}|{\Psi }_{\mathrm{el}}\left(t\right)⟩=\hat{H}|{\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}}⟩=\sum _p{C}_p^{\left(\lambda \right)}|{\varphi }_p⟩,$$

$$\varphi \left({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right)$$

$$\equiv$$

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

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

$$=$$

$$\displaystyle\Braket{\cdots abc\cdots\left missing}{\hat{F}}{\cdots abc\cdots}$$

$$\displaystyle\Braket{\cdots rbc\cdots\left missing}{\hat{F}}{\cdots abc\cdots}$$

$$\displaystyle\Braket{\cdots ars\cdots\left missing}{\hat{F}}{\cdots abc\cdots}=0.$$

$$|\varphi ⟩=|\cdots abc\cdots ⟩=-|\cdots acb\cdots ⟩=|\cdots cba\cdots ⟩.$$

$${\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{\phi_{p}\left missing}{\hat{F}}{\phi_{q}}=\sum_{A,B}^{N_{{A}}}\sum_{i_{A}=1}^{n_{{\rm AO}}{(A)}}\sum_{j_{B}=1}^{n_{{\rm AO}}{(B)}}D^{(p)}_{i_{A}}D^{(q)}_{j_{B}}\Braket{\chi_{i_{A}}\left missing}{\hat{F}}{\chi_{j_{B}}}$$

$$\hat{\rho }\left(\mathbf{r}\right)=\sum _k^N\delta \left(\mathbf{r}-{\mathbf{r}}_k\right)=\sum _k^N{\delta }_k(\mathbf{r}),$$

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

$$\hat{j}\left(\mathbf{r}\right)=\frac{1}{2}\sum _k^N\left({\delta }_k(\mathbf{r}){\hat{p}}_k+{\hat{p}}_k^{†}{\delta }_k(\mathbf{r})\right).$$

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

$$=$$

$$\int {\Phi }_{\lambda }^{\mathrm{CI}}\left({\mathbf{r}}^N\right)\hat{\rho }\left(\mathbf{r}\right){\Phi }_{\nu }^{\mathrm{CI}}\left({\mathbf{r}}^N\right)\mathrm{d}{\mathbf{r}}^N$$

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

$$=$$

$${\mathbf{j}}_{\lambda \nu }\left(\mathbf{r},t\right)=-\frac{i}{2}\sum _k(\int {\delta }_k(\mathbf{r})\left({\Phi }_{\lambda }^{\mathrm{CI}}\left({\mathbf{r}}^N\right){\vec{\nabla }}_k{\Phi }_{\nu }^{\mathrm{CI}}\left({\mathbf{r}}^N\right)\right)\mathrm{d}{\mathbf{r}}^N$$

$$-\int {\delta }_k(\mathbf{r})\left({\Phi }_{\nu }^{\mathrm{CI}}\left({\mathbf{r}}^N\right){\vec{\nabla }}_k{\Phi }_{\lambda }^{\mathrm{CI}}\left({\mathbf{r}}^N\right)\right)\mathrm{d}{\mathbf{r}}^N).$$

$$\sum _k\int {\delta }_k(\mathbf{r})\left({\Phi }_{\lambda }^{\mathrm{CI}}\left({\mathbf{r}}^N\right){\vec{\nabla }}_k{\Phi }_{\nu }^{\mathrm{CI}}\left({\mathbf{r}}^N\right)\right)\mathrm{d}{\mathbf{r}}^N=\sum _{p\ne q}{C}_p^{\left(\lambda \right)}{C}_q^{\left(\nu \right)}\left({\phi }_{\overline{r}}(\mathbf{r})\vec{\nabla }{\phi }_{\overline{a}}(\mathbf{r})\right).$$

$$\begin{array}{rl}\mathbf{y}\left(\mathbf{r},t\right)& ={\int }_0^t\mathrm{d}{t}^{\prime }\frac{\partial \rho \left(\mathbf{r},t^{\prime }\right)}{\partial t^{\prime }}=\rho \left(\mathbf{r},t\right)-\rho \left(\mathbf{r},0\right)\end{array}$$

$${\hat{\mu }}_r=-\int \mathbf{r}\hat{\rho }\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}.$$

$${\hat{\mu }}_v=-\int \hat{ȷ}\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}.$$

$${⟨{\left({\hat{\mu }}_v\right)}_r⟩}_{\lambda \nu }$$

$$|{\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),$$

$$|{\Psi }_{\mathrm{e}}⟩=\frac{1}{\sqrt{2}}\left(|{\Psi }_x⟩+i|{\Psi }_y⟩\right).$$

Full text

An Open-Source Framework for Analyzing $N$-Electron Dynamics: I. Multi-Determinantal Wave Functions

Vincent Pohl
, Gunter Hermann22footnotemark: 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.

Abstract

The aim of the present contribution is to provide a framework for analyzing and visualizing the correlated many-electron dynamics of molecular systems, where an explicitly time-dependent electronic wave packet is represented as a linear combination of $N$-electron wave functions. The central quantity of interest is the electronic flux density, which contains all information about the transient electronic density, the associated phase, and their temporal evolution. It is computed from the associated one-electron operator by reducing the multi-determinantal, many-electron wave packet using the Slater-Condon rules. Here, we introduce a general tool for post-processing multi-determinant configuration-interaction wave functions obtained at various levels of theory. It is tailored to extract directly the data from the output of standard quantum chemistry packages using atom-centered Gaussian-type basis functions. The procedure is implemented in the open-source Python program detCI@ORBKIT, which shares and builds upon the modular design of our recently published post-processing toolbox [J. Comput. Chem. 37 (2016) 1511]. The new procedure is applied to ultrafast charge migration processes in different molecular systems, demonstrating its broad applicability. Convergence of the $N$-electron dynamics with respect to the electronic structure theory level and basis set size is investigated. This provides an assessment of the robustness of qualitative and quantitative statements that can be made concerning dynamical features observed in charge migration simulations.

1 Introduction

The ultrafast evolution of transient electronic densities plays a central role in understanding the chemical reactivity and in predicting spectroscopic properties of molecules. With the recent advances in attosecond laser technologies, it has now become possible to indirectly observe the dynamics of electrons on their natural timescale [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Whereas direct observation of the electron flow remains elusive, its experimental reconstruction yields a wealth of information about charge migration in molecules, opening a wide range of new applications.

The theoretical description of electron dynamics has also greatly progressed over the last decades [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]. In particular, Time-Dependent Density Functional Theory (TDDFT) [48] holds lot of promises due to its computational efficiency and intuitive interpretation. Approaches based on $N$-electron wave functions such as Multi-Configuration Time-Dependent Hartree-Fock [23, 24, 28] or time-dependent Configuration Interaction (CI) [21, 49, 50, 29, 51, 39, 40, 52, 53, 54] offer an attractive alternative to density-based schemes. They share the common philosophy of representing an $N$-electron wave packet as a linear combination of spin-symmetrized Slater-determinants, which are constructed by exciting electrons from a single reference determinant. As such, the time-evolving wave packet is thus a multi-determinantal wave function. These methods are systematically improvable and converge to the same Full CI limit. Their major limitation is the high associated numerical cost, but this problem is mitigated by the ever increasing computational resources at our disposal.

The choice of an $N$-electron determinant basis for the electron dynamics has the important advantage of ensuring the $N$-representability of the wave packet at all times. It also reveals information about the dynamical build up of correlation in the transient electronic wave packet and its physical origin (particle-hole ($p$–$h$), two-particle–two-hole (2$p$–2$h$), $\dots$, excitations). On the other hand, it has the associated disadvantage of rendering interpretation of the electron dynamics less intuitive, which can be circumvented by using a proper set of visualization tools. Apart from the transient one-electron density, these need to include information about the evolution of the phase of the wave packet, which strongly depends on the electronic correlations. This complementary information is encoded in the electronic flux density, equivalent to the current density, from which a qualitative picture of the electron flow emerges naturally. The electronic flux density is a vector field that allows, at first glance, for a microscopic understanding of the mechanisms at work during charge migration processes.

In this contribution, we introduce a framework for analyzing and visualizing the correlated many-electron dynamics of molecular systems based on the reconstruction of the electronic flux density from a general multi-determinantal wave function. This requires a number of fundamental one-electron quantities, such as difference electronic densities, transient electronic flux densities, and transition dipole moments, that are not directly accessible from the output of standard quantum chemistry packages. Our initiative embraces the open-source molecular modelling philosophy [55] and builds upon the modular structure of our recently published quantum chemistry toolbox ORBKIT [56]. The latter is capable of computing a multitude of static electronic properties based on the data of electronic structure calculations from single-determinant wave function approaches, such as Hartree-Fock (HF) or Density Functional Theory (DFT) methods. Here, we extend the capabilities of ORBKIT to multi-determinant wave functions by exploiting its highly modular and easily comprehensible Python architecture. The new customized post-processing program, detCI@ORBKIT, can extract the data of multi-determinantal wave functions from various quantum chemistry programs using atom-centered Gaussian-type basis functions, and evaluate matrix elements of one-electron operators in the basis of $N$-electron eigenstates. The time-dependent quantities required for analyzing and visualizing the $N$-electron wave packet dynamics by means of the flux density are then calculated as linear combinations of the static matrix elements. This new tool will prove valuable to investigate a great number of charge migration processes. The present contribution also explores parameters that affect the qualitative and quantitative analysis of the electronic flux density, as applied to the electron dynamics in H${}_{3}^{+}$ and LiH.

The paper is organized as follows: Section 2 briefly describes the time-dependent configuration-interaction methodology and introduces general computational rules for computing one-electron matrix elements. The influence of the basis set size on the flux density is benchmarked in subsection 3.1 for the H${}_{3}^{+}$ test system. Subsection 3.2 investigates the influence of the electronic structure method on the qualitative features of the electron migration process. Concluding remarks are presented in Section 4. Atomic units are used throughout the paper ($\hbar=m_{e}=e=4\pi\varepsilon_{0}=1$), unless stated otherwise.

2 Computational Procedure and Theory

2.1 Time-Dependent Configuration Interaction

The electron dynamics of a molecular system can be described by solving the non-relativistic time-dependent electronic Schrödinger equation[57]

[TABLE]

The field-free Hamiltonian $\hat{H}$ for a system consisting of $N$ electrons and $N_{A}$ nuclei is written in the clamped nuclei approximation as

[TABLE]

Here, $r_{ij}$ is the distance between electrons $i$ and $j$, $Z_{A}$ is the charge number of nucleus $A$, and $r_{Ai}$ is the distance between nucleus $A$ and electron $i$. In general, the multi-particle time-dependent electronic wave function $\Ket{\Psi_{{\rm el}}\left(t\right)}$ can be formulated as a linear superposition of stationary electronic states $\Ket{\Phi_{\lambda}}$

[TABLE]

with $B_{\lambda}\left(t\right)$ as the time-dependent expansion coefficients of state $\lambda$. From a dynamical perspective, it is convenient to choose a basis of $N$-electron states that diagonalizes the field-free Hamiltonian at a given level of theory. This is the approach followed in the present paper.

Generally speaking, the time-independent $N$-electron wave function $\Ket{\Phi_{\lambda}}$ can be expressed in the terms of a configuration interaction (CI) expansion. CI methodologies are conceptually similar to other high-level wave function based methods, as will be discussed below. To ensure a proper description of electron correlation, the wave function is expanded in a complete set of configuration state functions[58]

[TABLE]

where the expansion parameters (or CI-coefficients) $C_{p}^{\left(\lambda\right)}$ are optimized variationally. In the present implementation, the orthonormal configurations $\Ket{\phi_{p}}$ are chosen as Slater determinants. These are defined as antisymmetrized products of one-electron spin orbitals $\Ket{\varphi_{a}}$

[TABLE]

Here, the Slater determinant is represented as a function of the spin and spatial coordinates $\mathbf{x}$ of the $N$ electrons, and $\left\{\varphi_{a}\varphi_{b}...\varphi_{c}\right\}$ are the occupied orthonormal molecular spin orbitals. In the CI-approach, the various Slater determinants are constructed by excitations of spin orbitals from a single reference state $\Ket{\phi_{0}}$. For example, exciting an electron from an occupied spin orbital $a$ to a virtual spin orbital $r$ from the reference state forms a singly excited determinant $\Ket{\phi_{a}^{r}}$. Accordingly, the CI wave function can be reformulated in terms of these excited determinants

[TABLE]

where $\left\{a,b,c\right\}$ are the occupied spin orbitals, and $\left\{r,s,t\right\}$ denote the virtual spin orbitals (i.e., unoccupied in the reference determinant). The exact ansatz including all possible excitations is referred to as Full CI approach. Since the number of conceivable excitations increases factorially with the number of electrons and orbitals, this is computationally very expensive and only feasible for small systems. To circumvent this bottleneck, two main strategies are pursued to reduce the Full CI expansion (cf. Eq. (11)): first, the truncation to a certain maximum rank of excitations, e.g., CI Singles (CIS) or CI Doubles (CID), and second, the restriction of the active space to a certain number of electrons in a specified number of orbitals, e.g., Multi-Configuration Self Consistent Field (MCSCF)[59]. In the latter scheme, the orbitals themselves appearing in Eq. (10) are variationally optimized in addition to the CI-coefficients. This yields a better representation of the correlation in a reduced orbital space, usually brought to the Full CI limit in an active space chosen close the HOMO-LUMO gap of the molecule. Quite importantly, all information required to reconstruct the orbitals and the $N$-electron eigenfunctions at a chosen level of theory are accessible from the output of standard quantum chemistry packages. Indeed, this information is used by many post-processing programs for computing orbital-derived quantities. We will now investigate how it is possible to use this knowledge to compute time-dependent wave function derived properties.

2.2 General Considerations on Expectation Values

Exploiting the structure of the time-dependent multi-determinant wave functions, Eqs. (3) and (4), the expectation value of any one-electron operator $\hat{F}$ can be expressed as

[TABLE]

In order to obtain the expectation value of the operator $\hat{F}$, one has to evaluate the matrix elements \Braket{\phi_{p}\left missing}{\hat{F}}{\phi_{q}} between two determinants (cf. Eq. (14)). For that purpose, the Slater-Condon rules[60, 61, 62], allow to express the respective matrix elements in terms of one-electron integrals in the spin orbital space. These rules can be summarized as follows \Braket{\phi_{p}\left missing}{\hat{F}}{\phi_{q}} for three general types:

Identical Slater determinants

[TABLE] 2. 2.

Two Slater determinants differing by a single spin orbital

[TABLE] 3. 3.

Two determinants which differ by two or more spin orbital

[TABLE]

A prerequisite for applying the Slater-Condon rules is the maximum coincidence principle, i.e., all common spin orbitals of both configurations appear at the same positions in the respective Slater determinant. This is achieved by permutation of the spin orbitals in one of the determinants, i.e., by interchanging columns in Eq. (9). This necessary re-ordering can change the sign due to the antisymmetric properties of determinants

[TABLE]

The resulting phase factor, $\left(-1\right)^{N_{{\rm P}}}$, can be determined by counting the required number of column permutations $N_{{\rm P}}$ to reach maximum coincidence. It must be accounted for when applying the Slater-Condon rules to compute the matrix elements in Eq. (14).

The computation of expectation values of any one-electron operator from a configuration interaction wave packet of the type Eq. (4) boils down to evaluate transition moments between spin orbitals. In computational chemistry, the spin orbitals are usually transformed to spin-free representations by integrating over the spin coordinates. We specialize here to the case, where these spatial molecular orbitals (MO) are defined in the framework of the MO-LCAO (Molecular Orbital - Linear Combination of Atomic Orbitals) ansatz. Specifically, an MO $\varphi_{a}$ is expanded using a finite set of atom-centered basis functions

[TABLE]

with $D^{(a)}_{i_{A}}$ as the $i_{A}$th expansion coefficient for MO $a$. The basis functions $\chi_{i_{A}}$ are atomic orbitals, 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)}$ denotes the number of atomic orbitals on atom $A$, with $N_{\rm AO}=\sum_{A=1}^{N_{A}}n_{{\rm AO}}{(A)}$. In the MO-LCAO representation, the transition moments between spin orbitals take the form

[TABLE]

All required information to reconstruct molecular orbitals, i.e., the MO-LCAO coefficients $D^{(a)}_{iA}$ and the atom-centered basis functions, can be found in the output of standard quantum chemistry program packages. The integrals in the atomic orbital basis are computed analytically with our Python post-processing toolbox ORBKIT for a wide range of one-electron operators [56].

2.3 The Electronic Continuity Equation

A widespread quantity for the visual representation of electronic motions in molecular systems is the time-dependent one-electron density $\rho\left(\mathbf{r},t\right)$.[63, 64, 65, 66] It remains a useful tool to characterize the correlated electron dynamics from multi-determinant wave functions of the form Eq. (4), and it can be computed as an expectation value of the associated one-electron density operator

[TABLE]

where $\delta\left(\mathbf{r}-\mathbf{r}_{k}\right)=\delta_{k}(\mathbf{r})$ is the Dirac delta distribution, $\mathbf{r}$ designates the grid of observation, and $\mathbf{r}_{k}$ refers to the position of electron $k$. The one-electron density admits a realistic interpretation as a probability fluid, which must satisfy a continuity equation of the form

[TABLE]

The vector field $\mathbf{j}\left(\mathbf{r},t\right)$ is the electronic flux density, often called current density, which contains information about electronic phase driving the spatial flow of the electron density. The associated operator can be written as

[TABLE]

Here, $\hat{p}_{k}=-i\vec{\nabla}_{k}$ is the momentum operator of an electron $k$, where $\vec{\nabla}_{k}$ is the gradient operator.

Using the time-dependent wave function of Eq. (3), the expectation value for the electron density reads

[TABLE]

with $\Psi_{{\rm el}}$ as a function of the spatial coordinates $\mathbf{r}^{N}$ of $N$ electrons. In Eq. (25), the matrix expression can be simplified by applying the Slater-Condon rules

[TABLE]

with $n_{\overline{a}}$ as the occupation number of MO $\overline{a}$. The over-line “–” denotes a formal excitation from the MO $\overline{a}$ to MO $\overline{r}$ taking the configuration state function of state $\nu$ as a reference. The corresponding expression for the electron flux density can be formulated as

[TABLE]

where $\mathbf{j}_{\lambda\nu}\left(\mathbf{r},t\right)$ is the transition electronic flux density from state $\lambda$ to state $\nu$

[TABLE]

Since the electronic states are real-valued, the diagonal terms $\mathbf{j}_{\lambda\lambda}\left(\mathbf{r},t\right)$, i.e., the adiabatic flux density [67, 68, 69, 70, 71], vanish. The same argument holds for all matrix elements in Eq. (28) involving identical determinants. Using the Slater-Condon rules, the integrals in Eq. (29) simplify to one-electron integrals over the spin orbitals

[TABLE]

The derivatives of the molecular orbitals are computed analytically using functions from the Python toolbox ORBKIT, with which both, the density and the electronic flux density, can then be projected on an arbitrary grid.

From the one-electron density, it is possible to derive another potentially useful quantity for the analysis of $N$-electron dynamics. The difference density $\mathbf{y}\left(\mathbf{r},t\right)$ describes the variation of the electron density within the time interval $\left[0,t\right]$ from a chosen initial condition[72, 73]. It is defined as the integral over time of the electron flow

[TABLE]

The difference density determines the number of electrons that have moved in and out of a specific volume element during a given laps of time. As such, it yields quantitative information that is complementary to both the electronic flux density and the electron flow.

The convergence of the continuity equation can be a priori estimated by inspecting closely related static quantities that relates expectation values of $\hat{\rho}\left(\mathbf{r}\right)$ and $\hat{\jmath}\left(\mathbf{r}\right)$. A potentially useful tool in this endeavor is the comparison of the dipole operator in length and in velocity gauge.[63] The former derives from the one-electron density and takes the form

[TABLE]

The latter is defined from the electronic flux density as

[TABLE]

As the transition moment between a given pair of states $\{\Phi_{\lambda},\Phi_{\nu}\}$, both can be directly related to each other via[74]

[TABLE]

The last expression can be used to estimate the quality of the level of theory and of the underlying basis of a given quantum chemical calculation. For a converged calculation, the transition moment $\Braket{\hat{\mu_{r}}}_{\lambda\nu}$ (Eq. (32)) and $\Braket{\left(\hat{\mu}_{v}\right)_{r}}_{\lambda\nu}$ (Eqs. (33) and (34)) must become identical for every pair of states involved in the dynamics.

2.4 Implementation Details

The computational rules described in the previous subsections are implemented in our new open-source Python toolbox detCI@ORBKIT, which is freely available at https://github.com/orbkit/orbkit/. The program requires a preliminary quantum chemistry calculation at a desired level of wave function theory, and using Gaussian-type atom-centered orbitals. Starting from the data of a determinantal CI-calculation, the program builds a library of transition moments and expectation values of one-electron operators projected on an arbitrary grid, to be used for analyzing the $N$-electron dynamics. The electronic ground state and all excited states serving as a basis for the subsequent dynamics must be computed at the same level of theory and using the same atomic basis. The toolbox detCI@ORBKIT is written in Python, which offers a broad set of efficient standard libraries and simplifies its portability on different platforms. The structure of the program can be summarized as follows:

A parser routine extracts the information about the molecular geometry, the atomic basis, the coefficients of the molecular orbitals and their energies, the MO occupation patterns for all desired electronic states, and the $N$-electron eigenstate energies. Currently, the program supports the MOLPRO, PSI4, GAMESS, and TURBOMOLE formats. A full and updated list can be found in the program documentation. 2. 2.

The molecular orbitals and the derivatives thereof are reconstructed from the atomic orbitals and projected on an arbitrary grid using the functionalities of ORBKIT. All integrals required for computing the matrix elements of the one-electron operators in the molecular orbital basis are computed analytically via the underlying atomic basis. These are stored in a Python list for later use. 3. 3.

A library of transition moments is built for each pair of $N$-electron stationary wave functions used in the time-dependent wave packet expansion (cf. Eq. (4)). The evaluation of transition moments of one-electron operator between two multi-determinant states, \Braket{\Phi_{\lambda}^{{\rm CI}}\left missing}{\hat{F}}{\Phi_{\nu}^{{\rm CI}}}, is performed in three steps:

1. (a)

The Slater determinants are compared to each other to determine all matrix elements \Braket{\phi_{p}\left missing}{\hat{F}}{\phi_{q}} involved. The configurations, where both CI-coefficients $\left\{C_{p}^{\left(\lambda\right)},C_{q}^{\left(\nu\right)}\right\}$ are larger than a user-defined threshold, are brought to maximum coincidence. The necessary number of orbital permutations determines the phase factor for the rearrangement. Since the number of Slater-determinants in a stationary state can become prohibitively high, the comparison and ordering routine is implemented in Cython[75] and can be executed on multiple processors. 2. (b)

From the occupation pattern two cases are identified: two identical determinants, and two determinants that differ by a single spin orbital. These build an integer list of important contributions. All one-electron integrals between two spin orbitals, \Braket{\varphi_{a}\left missing}{\hat{F}}{\varphi_{a}} and \Braket{\varphi_{a}\left missing}{\hat{F}}{\varphi_{r}}, which are necessary to calculate the expectation value of the operator, are loaded from the MO Python list generated by ORBKIT in step 2. 3. (c)

The transition moments \Braket{\Phi_{\lambda}^{{\rm CI}}\left missing}{\hat{F}}{\Phi_{\nu}^{{\rm CI}}} in Eq. (13) are calculated by weighting the integrals over the spin orbitals with the associated CI-coefficients $\left\{C_{p}^{\left(\lambda\right)},C_{q}^{\left(\nu\right)}\right\}$. In general, weighting the CI-coefficients and the phase factor with the grid-dependent molecular orbitals is the bottleneck of our routine, and significant effort was made on the optimization and parallelization of the associated Cython[75] modules. 4. 4.

The library can then be kept in memory or stored on disk for later use. The Python architecture allows direct visualization of the one-electron quantities using packages such as matplotlib[76] or Mayavi[77]. In addition, various file formats (e.g., standard Gaussian cube-files, hierarchal HDF5 data files, etc.) are supported, which enables visualization with external graphical programs, such as, VMD[78] or ZIBAmira[79].

The elements in the CI-determinant basis yield a matrix representation of the electronic Hamiltonian, which can be used to investigate the $N$-electron dynamics of a wave packet of the form Eq. (3). The analysis is performed by weighting with the time-dependent coefficients $B_{\lambda}{\left(t\right)}$ the expectation values and transition moments involved in a desired operator. These are the position representation for the one-electron density, Eq. (24), and the transition moments, Eq. (29), for the electronic flux density, Eq. (2.3). The dynamical program is not part of the standard detCI@ORBKIT implementation.

3 Numerical Examples

In this section, the capabilities of our toolbox to study the correlated electron dynamics in real-time are illustrated for different molecular systems. Two test molecules are studied to reveal the dependence of the computational procedure towards the quality of the electronic structure theory and influence of the basis set: the trihydrogen cation H${}_{3}^{+}$ and the lithium hydride molecule, LiH. A Python execution code for each of these examples and the data from the associated quantum chemical calculation are available in the program package. Note that a development version of detCI@ORBKIT was already used for several applications in the literature, see Refs. [80, 81, 82].

3.1 Basis Set Dependence of the Continuity Equation

We advocate using the time-dependent electron density $\rho\left(\mathbf{r},t\right)$ and time-dependent electronic flux density $\mathbf{j}\left(\mathbf{r},t\right)$ as complementary quantities for the analysis of correlated electron dynamics in molecular systems. The dependence of an $N$-electron dynamics on the underlying atomic basis set is an important convergence parameter that determines the quality of this analysis. To assess the robustness of the predictions concerning the time-dependent electronic flux density towards the basis set size, the trihydrogen cation H${}_{3}^{+}$ is studied using the minimal basis set, STO-3G, as well as a systematic series of Dunning-type basis sets, cc-pVXZ and aug-cc-pVXZ with X=D, T, Q, 5.[83, 84] In these calculations, the H${}_{3}^{+}$ molecule is chosen to retain the equilateral triangular equilibrium geometry ($r_{\rm HH}=0.87\,{\rm\AA}$) of the ground state $1\,{}^{1}A_{1}^{\prime}$ aligned in the $xy$-plane. To avoid artifacts coming from the electronic structure method, all calculations are performed at the Full CI level, using the open-source package PSI4[85].

The reasons for selecting H${}_{3}^{+}$ as a test system is threefold: First, the quantum chemical calculations can be performed at the Full CI limit for a large selection of basis sets due to the small number of electrons in H${}_{3}^{+}$. Second, its electronic structure is already well-studied due to its important role in interstellar chemistry.[86, 87, 88] This enables to compare with high-quality reference calculations.[89] Finally, initial conditions can be chosen such as to drive an interesting charge migration process that induces a periodic unidirectional circular current in the H${}_{3}^{+}$ plane. As was previously shown for other high-symmetry ring-shaped molecules [90, 82, 91], this can be achieved by a carefully chosen superposition of the ground state $1\,{}^{1}A_{1}^{\prime}$ and the degenerate state $1\,{}^{1}E^{\prime}$. In the basis of Full CI eigenstates used to study the $N$-electron dynamics, the field-free evolution of the system is known analytically at all times. In the particular case presented here, the time evolution of a wave packet consisting of two superposition states takes the following form

[TABLE]

where $\Ket{\Psi_{\rm g}}$ is the stationary wave function of the ground state $1\,{}^{1}A^{\prime}$, and $\Ket{\Psi_{\rm e}}$ denotes the stationary wave function of the degenerate excited state $1\,{}^{1}E^{\prime}$. The latter is chosen as a complex-valued linear combination of $1\,{}^{1}E_{x}^{\prime}$ and $1\,{}^{1}E_{y}^{\prime}$ with the relative phase set to $\eta=0$

[TABLE]

Here, $\Ket{\Psi_{x}}$ and $\Ket{\Psi_{y}}$ refer to the wave functions of state $1\,{}^{1}E_{x}^{\prime}$ and $1\,{}^{1}E_{y}^{\prime}$, respectively. The axis labels correspond to the molecular orientation as given in the quantum chemistry program. It was shown that it is possible to prepare such a wave packet by electronic excitation of the ground state using a circular polarized laser field.[82] The time-dependent electron density associated with this wave packet takes the form

[TABLE]

with $\Delta E=E_{\rm{e}}-E_{\rm{g}}$. The first term on the right-hand-side describes a static contribution to the one-electron density, where $\left\{\rho_{\rm g},\rho_{x},\rho_{y}\right\}$ are the respective contributions from the ground state $1\,{}^{1}A_{1}^{\prime}$ and the excited states $1\,{}^{1}E_{x}^{\prime}$ or $1\,{}^{1}E_{y}^{\prime}$. The evolution of the electronic wave packet is driven by the transition electron densities, $\left\{\rho_{{\rm g}x},\rho_{{\rm g}y}\right\}$, which are obtained by resolving Eq. (24) in the basis of Full CI eigenstates. Similarly, the time-dependent electronic flux density obtained from Eq. (2.3) reads

[TABLE]

where $\left\{\mathbf{j}_{{\rm g}x},\mathbf{j}_{{\rm g}y},\mathbf{j}_{xy}\right\}$ stand for the transition electronic flux densities between the ground and the excited states. For a longer derivation of Eqs. (37), (38), the reader is referred to previous work. [82]

In order to determine the dependency of the electron density and the transition flux density towards the quality of the basis set, we make use of the continuity equation. The time derivative of the electron density $\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$ (left-hand side of Eq. (22)) will be hereafter referred to as “electron flow”, to differentiate this quantity from the divergence of the electronic flux density $-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$ (right-hand side of Eq. (22)). Both quantities represent expressions for the electron flow and should converge to the same value by definition. Since Full CI calculations are exact for a given basis set, any difference between both solely stems from the quality of the basis set. For clarity, only a few of the basis sets mentioned above are compared here. These include the minimal basis set STO-3G, and three correlation consistent basis sets, cc-pVTZ, aug-cc-pVTZ, and aug-cc-pV5Z. The results of Fig. 1 show the electron flow $\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$ (central panels), the transient electronic flux density $\mathbf{j}\left(\mathbf{r},t\right)$ (left panels), and its divergence $-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$ (right panels) for the superposition of the ground state $1\,{}^{1}A_{1}^{\prime}$ and the excited state $1\,{}^{1}E^{\prime}$ at $t=\tau/4$ ($\tau=h/\left(\Delta E\right)$). The red (blue) areas in the central and right panels represent regions of instantaneous decrease (increase) of the electron density. At first glance, it can be observed that the qualitative features of the electron flux density and the electron flow are very robust towards the basis set quality. That is, the electrons migrate from the lower left hydrogen atom to the lower right one. Further, the s-character of the three orbitals involved in the H${}_{3}^{+}$ bond can be quite easily recognized. This character is retained when using more complete basis sets for the analysis of the electron flow based on the density. As a complementary analysis tool, the electronic flux density maps reveal that the density migrates from one hydrogen to the other along along the bond at this time step. This somewhat counterintuitive feature is recognized for all basis sizes.

Nonetheless, several differences between $\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$ and $-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$ can be identified across the different basis sets. On the one hand, the electron flow and the vorticity in the electronic flux density obtained with the minimal basis set (STO-3G, top panels) are much larger than the results from the Dunning-type basis sets. This is due to the small number of basis functions, which overestimates the contribution of the $s$-orbitals to the total density. The electron flow appears to be converged already at the cc-pVTZ level (second line). When comparing $\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$ and $-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$, artifacts can be observed around the hydrogen nuclei. Since the electron density $\mathbf{\rho}\left(\mathbf{r},t\right)$ is rather insensitive towards the basis set quality, these artifacts merely occur when computing the divergence of the electronic flux density, $-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$. Using the cc-pVDZ basis (not shown), the nodal structures at the nuclei are largest and they disappear slowly as the number of basis functions increases. The reason for the slow convergence of the divergence of the electronic flux density in the present example is that it is dominated by the derivative of the density at the nuclei. Since the cusp at a nucleus is poorly represented using Gaussian-type atomic orbitals, the derivative of the transient flux density is mostly affected.

The same phenomenological robustness is observed for the time-evolution of the electronic flux density $\mathbf{j}\left(\mathbf{r},t\right)$ and the electron flow, $\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$. In Fig. 2, the time-dependent flux densities (arrows) are superimposed on the time-dependent electron flow, where the red (blue) areas indicate regions of decreasing (increasing) density. They are depicted at characteristic times during half the period $\tau$ of the charge migration process, $0\leq t\leq\tau/2$. The period $\tau$ is related to the energy difference $\Delta E$ as follows $\tau=h/\Delta E$. The results are shown for the cc-pVTZ (top panels) and the aug-cc-pV5Z basis set (bottom panels), for which the energy difference is $\Delta E=19.36\,{\rm eV}$ and $\Delta E=19.33\,{\rm eV}$, respectively. Both tools correctly predict a clockwise circular migration of the electron for the transition between state $1\,{}^{1}A_{1}^{\prime}$ and state $1\,{}^{1}E^{\prime}$. Neither qualitative nor quantitative differences between the two basis sets can be recognized. From Eq. (38), it can be observed that the electronic flux density has three components: a static ring current $\mathbf{j}_{xy}$, and two alternating polarized components $\mathbf{j}_{gx}$ and $\mathbf{j}_{gy}$. The latter two create the asymmetric pattern observed in the flux density, which shows that the density hops from atom to atom along an inward curved path (see, e.g. at $t=\tau/8$). Quite importantly, this information cannot be obtained from neither the density nor the electron flow.

In order to accurately determine the influence of the basis set size on the electron flow ($\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$, solid lines) and the divergence of the electronic flux density ($-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$, dashed lines), both quantities are illustrated in Fig. 3 as a function of the polar angle $\alpha$ at $t=\tau/4$. The functions are obtained by projecting the electronic density on a cylindrical grid, $\{r,\alpha,z\}$, and integrating both sides of the continuity equation over $r\,{\rm d}r\,{\rm d}z$. The grid representation is here again produced using ORBKIT and the integrals are performed via cubature [92, 93, 94]. The vertical dashed black lines denote the position of the nuclei. As previously observed, the minimal basis STO-3G is seen to poorly satisfy the continuity equation, as the discrepancies between the two quantities (gray curves) remain large over the whole domain $[0,2\pi]$. In comparison, the Dunning-type basis sets perform very well over the whole angular range already at the cc-pVDZ level. A look at the angular electron flow in the vicinity of the nuclei (cf. inset in Fig. 3) reveals that the plotted quantities ($\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$ and $-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$) – and with them the electronic continuity equation Eq. (22) – are quantitatively converged from the cc-pVTZ basis set level.

A further quantitative measure for the basis set convergence can be obtained by comparing the dipole moment in length gauge (cf. Eq. (32)) with the corresponding one calculated from the dipole moment in velocity gauge (cf. Eq. (34)). Both quantities must be equal for a converged basis set in a Full CI calculation. In Tab. 1, the $x$-components of the transition dipole moment for the two gauges are listed, along with the associated excitation energies for the $\Ket{\Psi_{g}}\to\Ket{\Psi_{x}}$ state transition and the total number of atomic basis functions $N_{\rm AO}$ for each basis set. The number of Full CI configurations for each state is obtained from the determinental CI output of the quantum chemistry program. This number is smaller for the ground state because all single excitations are projected out. The degenerate excited state exhibits a small splitting due to the use of the D2h abelian symmetry group in the calculation. Apart from the poor STO-3G basis, all results obtained from correlation consistent basis sets are quite homogeneous and yield a smooth convergence towards the literature benchmark value. As was the case for the electronic flux density, the cc-pVTZ values appear to be converged to a sufficient accuracy to allow for quantitative analysis. A similar level of accuracy on the energy and the dipole moment is obtained with a marginally smaller basis using diffuse functions, at the aug-cc-pVDZ level. Note that increasing the basis size also increases the number of Full CI configurations (see last columns of Tab. 1), which improves the variational description of the molecular orbitals and the $N$-electron wave functions simultaneously. These observations are complemented by Fig. 4, which shows the dipole moments and excitation energies as a function of the basis set size, ordered according to the number of functions. Here, the excitation energies are plotted as the difference energy to a highly accurate reference energy. The minimal basis set is excluded for clarity. From the inset in Fig. 4, it can be deduced that diffuse functions play a more important role for the energy convergence than adding basis functions with higher angular momenta. Interestingly, this does not apply to the convergence of the two dipole moments in length gauge, Eqs. (32) and (34), which converge monotonically with the number of basis functions.

3.2 Impact of the Electronic Structure Method

As a second example, the impact of different electronic structure methods on the time-dependent electron density and electronic flux density is investigated. In particular, both quantities are computed on the basis of a Full CI calculation and compared with those obtained from a CI Singles calculation, from a restricted active space configuration interaction (RASCI), and from Complete-Active Space Self-Consistent Field (CASSCF) calculations. The details of the active space are given below. These methods build a hierarchy of electronic structure theories, where the description of electron correlation is improved systematically by considering either a larger active space or a higher degree of excitations. Since a Full CI calculation is used as a reference, a well-studied four-electron molecule, the heteronuclear polar lithium hydride LiH, is chosen as a test system.[95, 96, 97, 98, 99, 100, 101, 102] We are particularly interested in the charge transfer state A${}^{1}\Sigma^{+}$, which is optically accessible from the electronic ground state X${}^{1}\Sigma^{+}$ and lies in the Franck-Condon region. By using a so-called $\pi/2$-pulse, it is possible to create a superposition state as in Eq. (35), where $\Ket{\Psi_{\rm g}}$ is the stationary wave function of the ground state and $\Ket{\Psi_{\rm e}}$ that of the charge transfer state.

The one-electron density associated with this superposition state evolves in time according to

[TABLE]

where $\Delta E=E_{{\rm e}}-E_{{\rm g}}$. The static one-electron density of the ground state X${}^{1}\Sigma^{+}$ ($\rho_{\rm g}\left(\mathbf{r}\right)$) and charge transfer state A${}^{1}\Sigma^{+}$ ($\rho_{\rm e}\left(\mathbf{r}\right)$), as well as the static transition density ($\rho_{\rm ge}\left(\mathbf{r}\right)$), are computed using Eq. (26). Similarly, the electronic flux density takes the simplified form

[TABLE]

where the transition electronic flux density $\mathbf{j}_{\rm ge}\left(\mathbf{r}\right)$ is computed using Eq. (30). Note that, contrary to the previous example, the electronic flux density does not have a time-independent current term (cf. Eq. (38)). To complement the analysis of electron migration, the difference density[72, 73] is calculated for the specific superposition state from Eqs. (31), (39), and the choice of phase $\eta=\pi$ as

[TABLE]

As the electron flow itself, it is found to be independent of the static densities of the ground and excited states.

For each electronic structure method, a single-point calculation at the ground-state equilibrium geometry of LiH ($r_{\rm Li-H}=1.63\,{\rm\AA}$) are performed with PSI4[85] using an aug-cc-pVTZ basis set. On this basis, the time-dependent electronic flux densities and density differences are calculated at representative time steps during one period $\tau$ ($\tau=h/\left(\Delta E\right)$) of the charge migration: $\tau/4$, $\tau/2$, and $3\tau/4$. This choice allows for a direct comparison of the dynamics despite the different transition energies, and hence the different timescales, found using the various methods. The resulting $\mathbf{j}\left(\mathbf{r},t\right)$ and $\mathbf{y}\left(\mathbf{r},t\right)$ are depicted in Fig. 5. The four horizontal panels show the results for the different levels of electronic structure methods: (a) CIS, (b) RASCI(2,5), (c) CASSCF(2,5), and (d) Full CI. The energy difference between the ground and charge transfer states are found to be $\Delta E_{\rm CIS}=4.04\,{\rm eV}$, $\Delta E_{\rm RAS(2,5)}=4.13\,{\rm eV}$, $\Delta E_{\rm CAS(2,5)}=3.41\,{\rm eV}$, and $\Delta E_{\rm FCI}=3.56\,{\rm eV}$. By definition, $\mathbf{j}\left(\mathbf{r},t\right)$ and $\mathbf{y}\left(\mathbf{r},t\right)$ are zero at $t=0$ and $t=\tau$ and therefore, not depicted here.

In the example presented, the Full CI method serves as a reference, since it provides the exact solution for the correlated electronic wave function within the limitation of a finite basis set. We will thus first proceed to a qualitative analysis of the associated electronic flux density and density difference to highlight their main characteristics. The regions of electron density depletion are denoted in red and the electronic flux density is depicted using streamlines, an alternate representation that connects the arrows of the vector field to give the impression of a fluid in motion. Phenomenologically, an electron migration from the hydrogen atom to the lithium atom is observed during the first half period $\tau/2$, and the reverse process for the second half (cf. Fig. 5). This qualitative observation is indeed expected for the superposition of the more ionic ground state X${}^{1}\Sigma^{+}$ with the first excited state A${}^{1}\Sigma^{+}$ which has a more covalent character. As can be seen from the density difference, the electron depletion region surrounding the hydrogen atom (red areas) are much more diffuse than the area of charge concentration (gray areas) on the lithium. Contrary to the electron flow (not shown), the density difference retains the same sign in this superposition state at all times. This implies that an electron migrates to the lithium and back, leaving a hole on the hydrogen atom and filling it subsequently. In the earlier stages of the propagation (left panels of Fig. 5), the electronic flux density exhibits a large vorticity around the lithium atom. The charge is transferred indirectly from the hydrogen to the lithium, where the electron density flow forms a torus and the density is enriched from behind. This finding is in agreement with previous theoretical studies on similar molecules [99, 98, 102]. The electronic flux density offers the main advantage of revealing all mechanistic features of the electron flow at first glance.

The CIS method [103] provides a computationally cheap alternative to Full CI that only contains single excitations from the reference wave function (cf. Eq. (11)). Despite some well-known limitations[104, 105], this simple approach is generally expected to correctly recover the qualitative character of molecular excited states.[106, 107] This can likewise be confirmed by comparing the results for the LiH molecule (cf. Fig. 5 upper panel) with the Full CI benchmark (cf. Fig. 5 lower panel). It is striking that all features of both the flux density and of the electronic difference density, can be captured at the CIS level. In particular, the toroidal structure of the time-dependent electronic flux density is very well captured. However, quantitative differences are discernible for the density difference (see, e.g., the central panels). In the present case, the main effect seems to be the localization of the positive charge closer to the hydrogen nucleus. These discrepancies can be attributed to an insufficient representation of electron correlation at the CIS ansatz.

To investigate this effect further, alternative determinant-based calculations were performed at the CASSCF level of theory[108, 109, 110, 111]. This special form of the MCSCF method, which is briefly explained in section 2, incorporates higher-order excitations in the description of the correlated wave function. While the Full CI and the CIS scheme are very straightforward to set up, the choice of the active space for a CASSCF calculation is an art in itself. To design a first active space, an educated guess can be obtained by considering the molecule in the dissociation limit. In the ground state at the dissociation limit, the lithium atom is found in the configuration 1s22s in the ${}^{2}S$ state and the hydrogen atom is in the ${}^{2}S$ state (1s1), which yield the following degenerate molecular states: a singlet ${}^{1}\Sigma^{+}$ and a triplet ${}^{3}\Sigma^{+}$ state. These states form the lowest covalent dissociation limit. The second lowest atomic excitation is the transition of the lithium atom to the 1s22p configuration in the ${}^{2}P^{o}$ state. This corresponds to four molecular states: a singlet ${}^{1}\Sigma^{+}$ state, a singlet ${}^{1}\Pi$ state, a triplet ${}^{3}\Sigma^{+}$ state, and a triplet ${}^{3}\Pi$ state, which represent the second lowest dissociation limit. In the dissociation limit, only the 2s and 2p orbitals of the lithium are required to describe difference between the ground state X${}^{1}\Sigma^{+}$ and the first singlet excited state A${}^{1}\Sigma^{+}$. Consequently, the minimal active space consists of two active electrons in five molecular orbitals with the 1s2 of Li as core orbitals. Note that this analysis does not strictly apply for the molecule at the ground-state equilibrium geometry, but it serves as an initial guess.

According to the prescription above, the LiH molecule is first calculated using the minimal active space while keeping the orbitals frozen. This restricted active space configuration interaction ansatz [112] is labeled according to the CASSCF same notation as RASCI(2,5). The resulting flux density and the associated electronic difference density are depicted in the second row of Fig. 5. Again, the results for the RASCI(2,5) calculation correctly recover the qualitative aspects of the electron redistribution process, in particular the large vorticity of the toroidal vector field surrounding the lithium atom (see left and right panels). Nonetheless, the region of density depletion in the electronic difference density is found to be much more localized at the hydrogen nucleus compared to Full CI. The region of density enrichment to the left of the lithium atom at $t=\tau/2$ has also a smaller spatial extent than in the benchmark. These are strong indications that electron correlation in LiH is poorly described using RASCI(2,5).

To improve the description of the static correlation, the molecule is recalculated at the state-averaged CASSCF(2,5) level of theory. In general, state-averaging is required to simultaneously calculate degenerate states of different symmetries on the basis of a single set of optimized molecular orbitals. In our example, this single set is prerequisite to apply the Slater-Condon rules. As can be seen from the third row of Fig. 5, this setup yields both a qualitative and quantitative agreement with the reference Full CI results but with a fraction of determinants necessary to describe the correlated wave function. Provided a proper active space can be constructed for a given molecule, this can amount to very significant computational savings. Even for this small test system, the Full CI approach comprises 9428 Slater determinants, while only 11 are required for the CASSCF(2,5) wave function ansatz of the first excited state A${}^{1}\Sigma^{+}$. This number also compares advantageously to the 80 determinants required at the CIS level. As a bottom line, it appears that qualitative features of the electron migration process can be obtained at even relatively crude levels of electronic structure theory (e.g., CIS), but a careful treatment of electron correlation is required for quantitative predictions.

4 Conclusions

In this paper, we have introduced a general framework for post-processing determinant-based configuration-interaction wave functions. The primary goal of this open-source project is to develop a tool for the characterization and analysis of correlated electron dynamics in molecular systems, where a wave packet is expanded using static $N$-electron wave functions. The procedure relies on the numerical determination of transition moments of a set of one-electron operators, which yields a time- and space-resolved picture of the $N$-electron dynamics. These include transition densities, the electronic flux density, and various derived observables. All quantities required to reconstruct the multi-determinant wave functions are extracted from the output of standard quantum chemistry packages using Gaussian-type atom-centered basis sets. The entire procedure is implemented in a novel Python program detCI@ORBKIT which extends the functionalities of the post-processing toolbox ORBKIT. The latter calculates molecular electronic properties from the data of single-determinant wave functions which is also extracted from quantum chemical calculations. The new procedure is constructed so that it can principally evaluate transition moments of any one-electron operator, by taking advantage of Slater-Condon rules to drastically reduce the numerical effort. Emphasis was put on the general applicability, the parallelization of computationally demanding steps, and the easy visualization of the results.

In the application examples, we have demonstrated that the selected set of one-electron quantities is suitable to characterize the correlated electron dynamics for molecular systems in real time. In particular, analysis of the electron flux density reveals microscopic details about the motion and flow of the electrons during the investigated dynamical processes at first glance. Its qualitative analysis has proven very robust towards the choice of electronic structure theory method and the quality of the underlying atomic basis set, with the exception of the minimal basis STO-3G. Comparison of the electron flow (the time derivative of the electron density, $\partial\mathbf{\rho}\left(\mathbf{r},t\right)/\partial t$) with the divergence of the electronic flux density $-\vec{\nabla}\cdot\mathbf{j}\left(\mathbf{r},t\right)$ reveals a slow convergence of the continuity equation, Eq. (22), with respect to the basis size. It thus appears preferable to base quantitative predictions on observables derived from the electron density, in particular on the electron flow and its integral over time, the electronic difference density. In this respect, it appears that an accurate description of electron correlation is of primal importance. The tools advocated here appear as complementary for the analysis of $N$-electron dynamics.

5 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. The funding of the Deutsche Forschungsgemeinschaft through the Emmy-Noether program (project TR1109/2-1) and of the Elsa-Neumann foundation of the Land Berlin are also acknowledged. Furthermore, we thank Jhon Pérez-Torres for fruitful discussions during the initiation of this project and Prof. Jörn Manz for encouraging this paper.

6 Keywords

Correlated Electron Dynamics, Slater-Condon Rules, Multi-Determinant Wave Function, Electronic Flux Density, Electron Density, Electronic Difference Density, Electronic Current Density

7 TOC

Unraveling correlated electron dynamics: We introduce an open-source Python framework to post-process determinant-based configuration-interaction data from standard quantum chemistry packages. The procedure builds a library of transition moments of selected one-electron operators. The library can be used to visualize and analyze the time-evolution of a molecular system, represented as a time-dependent linear combination of multi-determinantal wave functions.