An Open-Source Framework for N-Electron Dynamics: II. Hybrid Density Functional Theory/Configuration Interaction Methodology
AAn Open-Source Framework for N -Electron Dynamics: II. HybridDensity Functional Theory/Configuration Interaction Methodology Gunter Hermann ∗†‡ , Vincent Pohl †‡ , and Jean Christophe Tremblay † March 13, 2018
In this contribution, we extend our framework for analyzing and visualizing correlated many-electrondynamics to non-variational, highly scalable electronic structure method. Specifically, an explic-itly time-dependent electronic wave packet is written as a linear combination of N -electron wavefunctions at the configuration interaction singles (CIS) level, which are obtained from a referencetime-dependent density functional theory (TDDFT) calculation. The procedure is implemented inthe open-source Python program detCI@ORBKIT , which extends the capabilities of our recentlypublished post-processing toolbox [ J. Comput. Chem. (2016) 1511]. From the output of standardquantum chemistry packages using atom-centered Gaussian-type basis functions, the frameworkexploits the multi-determinental structure of the hybrid TDDFT/CIS wave packet to compute fun-damental one-electron quantities such as difference electronic densities, transient electronic fluxdensities, and transition dipole moments. The hybrid scheme is benchmarked against wave functiondata for the laser-driven state selective excitation in LiH. It is shown that all features of the elec-tron dynamics are in good quantitative agreement with the higher-level method provided a judiciouschoice of functional is made. Broadband excitation of a medium-sized organic chromophore furtherdemonstrates the scalability of the method. In addition, the time-dependent flux densities unravelthe mechanistic details of the simulated charge migration process at a glance. ∗ corresponding author † Institut f¨ur Chemie und Biochemie, Freie Universit¨at Berlin, Takustraße 3, 14195 Berlin, Germany ‡ These authors contributed equally to this work. a r X i v : . [ phy s i c s . c h e m - ph ] A p r Introduction
Unraveling the flow of electrons inside a molecule out of equilibrium is key to understand its reactiv-ity. Since the pioneering laser experiments by Zewail and co-workers [1, 2], the development of newlight sources has now granted access to the indirect observation of electron dynamics on its naturaltimescale. To shed light on the mechanistic details of this attosecond dynamics, accurate theoreticalmethods are required that capture the subtle details of the transient electronic structure evolution.Various approaches based on explicitly time-dependent density functional theory (TDDFT) andwave 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 forultrafast dynamics in strong laser fields. On the other hand, the advantages of wave function-basedmethods in terms of convergence become rapidly compensated by their unfavorable computationalcost. Further, the intuitive picture of electrons flowing on a molecular skeleton can become blurredby correlation effects between the N particles.This contribution is motivated by the need for a robust, scalable wave function method toinvestigate ultrafast N -electron dynamics in systems of large dimension. The method we advocateis based on a combination of linear-response time-dependent functional theory (LR-TDDFT) andconfiguration interaction singles (CIS) methodologies, as was introduced recently [3–5]. In principle,the method is similar to the well-established CIS ansatz, with the exception that the energies andthe pseudo-CIS eigenvectors are obtained from a reference LR-TDDFT calculation. This allowsto improve the energetic properties of the states while keeping a simple electron-particle pictureto describe the transient N -electron wave packet. This TDDFT/CIS hybrid formalism inheritsthe qualities of both underlying methods and ensures the N -representability of all reduced densitymatrices, at all times and under all laser conditions. The ensuing N -electron dynamics remainsmarred by the non-intuitive interpretation of quantities beyond the density itself.Recently, we demonstrated that correlated electron dynamics can be accurately described bymeans of the electronic flux density operator and derived one-electron properties. We introducedan open-source framework [6] to post-process multi-determinantal configuration interaction wavefunctions directly from the output of standard quantum chemistry packages. It thus becomes pos-sible to reconstruct the transient N -representative one-electron density and current density (fluxdensity) using a library of transition moments calculated from the multi-determinantal configura-tion interaction wave functions, yielding an intuitive tool for visualizing and analyzing the corre-lated 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 formalismto the TDDFT/CIS hybrid formalism mentioned above, which should retain the qualities of thewave 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 thedescription of the analysis toolset based on the flux density. The application section reports on2enchmark calculations on the LiH molecule and the demonstration of the scalability of the schemeby investigation of the broadband excitation in an organic chromophore. The findings are sum-marized in the conclusion section. Unless otherwise stated, atomic units are used throughout themanuscript (¯ h = m e = e = 4 πε = 1). The evolution of the electronic state of a molecular system obeys the time-dependent Schr¨odingerequation [7], which can be written in the clamped nuclei approximation i ∂∂t | Ψ el ( t ) (cid:105) = (cid:16) ˆ H − ˆ µ · (cid:126)F ( t ) (cid:17) | Ψ el ( t ) (cid:105) . (1)The interaction of the molecular dipole ˆ µ 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 Hamilto-nian reads ˆ H = − N (cid:88) i =1 ∇ i + N (cid:88) i =1 N (cid:88) j>i r ij − N (cid:88) i =1 N A (cid:88) A =1 Z A r Ai , (2)where r ij = | (cid:126)r i − (cid:126)r j | is inter-electronic Coulomb repulsion, and r Ai is the distance between the i thelectron and nucleus A of charge Z A . In this work, an electronic wave packet | Ψ el ( t ) (cid:105) satisfyingEq. (1) is expressed as a linear superposition of stationary electronic states | Φ λ (cid:105)| Ψ el ( t ) (cid:105) = (cid:88) λ B λ ( t ) | Φ λ (cid:105) . (3)Here, B λ ( t ) are the expansion coefficients of state λ , which describe the time-evolution of the wavepacket. For molecules in strong laser fields, a large number of stationary electronic states is requiredto offer a proper description of the N -electron dynamics. The equations of motion for the coefficientsin Eq. (1), associated with the basis set expansion Eq. (3), can be integrated numerically.In the time-dependent configuration interaction methodology, the stationary electronic statesare chosen as linear combinations of excited configuration state functions (cid:12)(cid:12) Φ CI λ (cid:11) = C ( λ )ref | φ ref (cid:105) + (cid:88) ar C r ( λ ) a | φ ra (cid:105) + (cid:88) abrs C rs ( λ ) ab | φ rsab (cid:105) + . . . . (4)The expansion parameters C ( λ ) are associated with the formal excitation of a reference configuration, | φ ref (cid:105) , from occupied orbitals { a, b, c } to virtual orbitals { r, s, t } . Including all possible excitationsleads to the exact Full CI limit. The reference and excited configurations are defined as Slaterdeterminants, which builds antisymmetrized products of one-electron spin orbitals | ϕ a (cid:105) . 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 at3 given level of theory. The matrix elements of the dipole operator can be computed from theknowledge of these eigenfunctions, which serve as a basis for the variational representation of themolecule-field interaction.For large molecules, it is customary to truncate the CI expansion to a chosen maximum rankof excitations (e.g., CI Singles or CI Singles Doubles) in order to reduce the number of possible ex-cited configurations. Unfortunately, this often compromises the energetic description of the excitedstates. To circumvent this limitation while keeping the problem computationally tractable, Sonkand Schlegel [3] first recognized that only excitation energies and transition dipole moments arerequired to perform TDCI simulations. These can be obtained from linear-response time-dependentdensity functional theory (LR-TDDFT). To generalize this approach, it was proposed to use thesolutions of the LR-TDDFT calculation to generate a basis of pseudo-CI eigenstates [4, 5]. Allrequired information for a TDCI simulation is thus available from the output of standard quantumchemistry programs, provided excitations are performed from the ground state.According to the Runge-Gross theorem, it is possible to recast the N -electron Schr¨odingerequation and calculate all observables from the sole knowledge of the one-electron density. Usingthe Kohn-Sham ansatz for the density, the N -electron time-dependent Schr¨odinger equation can bemapped onto a one-electron equation for the orbitals i ∂ϕ a ( r , t ) dt = (cid:18) − ∇ v KS ( r , t ) (cid:19) ϕ a ( r , t ) . (5)The time-dependent Kohn-Sham potential v KS ( r , t ) contains the classical electrostatic interaction( v Hartree ( r , t )), an external potential ( v ext ( r , t )), and an exchange-correlation contribution ( v xc ( r , t )),i.e., v KS ( r , t ) = v Hartree ( r , t ) + v ext ( r , t ) + v xc ( r , t ) . (6)In explicit TDDFT, the Kohn-Sham potential is usually assumed to be local in time. A celebratedsuccess of TDDFT comes from its linear-response formulation, which allows to accurately computespectral properties of large molecules. For this endeavor, the response kernel of the electron densityto an external weak, long wavelength perturbation can be evaluated from the electric susceptibilityof the ground state. The search for the poles of the response function can be recast as an eigenvalueproblem of the form (cid:34)(cid:32) A BB † A † (cid:33) − ω (cid:32) − I I (cid:33)(cid:35) (cid:32) XY (cid:33) = − (cid:32) δ v δ v † (cid:33) , (7)where δ v is the response of the system state to the perturbation, δv ext ( r , t ). The elements of matrices A and B are obtained from the orbital energies and integrals over the exchange-correlation kernel,see Eq. (6). At the resonance frequencies ω , where the response vanishes ( δ v = 0), the solution ofthe Casida Eq. (7) yields simultaneously the excitations and de-excitations amplitudes, X and Y .In the present work, we make use of the fact that these are usually given in the output of standardquantum chemistry programs, together with the excitation energies and the oscillator strengths.4rom Eq. (7), it is possible to define pseudo-CI Singles eigenvectors in the Tamm-Dancoff ap-proximation, which consists in neglecting the off-diagonal blocks B . This procedure can alter thequality of the energetic properties of the excited states. On the other hand, the dominant charac-ters present in the pseudo-CI eigenvectors are often not strongly affected by this approximation. Inthe TDDFT/CI procedure, we thus advocate using directly the transition energies and amplitudesobtained from a LR-TDDFT calculation to take advantage of the full solution of Eq. (7) and toobtain a good energetic description of the excited states. A separate Tamm-Dancoff calculationmay be used to confirm the character of the excited states. The LR-TDDFT excitation amplitudesare then re-orthonormalized using a modified Gram-Schmidt procedure to define a pseudo-CI basisfor the TDCIS dynamics. All properties not directly deriving from the energies can be subsequentlycalculated at the CIS level of theory using the orbitals and the pseudo-CI eigenstates, which aretreated as configuration interaction singles expansions. Note that the Slater determinants { φ ref , φ ra } are constructed from Kohn-Sham orbitals. Importantly, all the information required to reconstructthese KS-orbitals and the N -electron pseudo-eigenfunctions are directly accessible from the outputof standard quantum chemistry packages. As a consequence, only the evaluation of one-electronintegrals is required to generate a library of molecular properties and transition moments of variousone-electron operators, which can be used to characterize the properties of transient wave packets,as explained below. For the analysis of the N -electron dynamics, we propose using a set of tools composed from theone-electron density, ρ ( r , t ), and the associated electronic flux density, j ( r , t ). These are related bythe electronic continuity equation ∂∂t ρ ( r , t ) = − (cid:126) ∇ · j ( r , t ) . (8)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. Theone-electron density can be used to define the electron flow, ∂∂t ρ ( r , t ), as the left-hand-side of thecontinuity equation. The difference density, y ( r , t ), is a widespread quantity used for visualizationpurposes, and it can be obtained by integrating the electron flow from a chosen initial condition ρ ( r , y ( r , t ) = (cid:82) t d t (cid:48) ∂ρ ( r ,t (cid:48) ) ∂t (cid:48) = ρ ( r , t ) − ρ ( r , . (9)We will resort to both quantities in later analyses.In operator form, the one-electron density and the electronic flux density respectively readˆ ρ ( r ) = (cid:80) Nk δ ( r − r k ) = (cid:80) Nk δ k ( r ) (10)ˆ j ( r ) = (cid:80) Nk (cid:16) δ k ( r )ˆ p k + ˆ p † k δ k ( r ) (cid:17) , (11)5here r is an observation point, δ ( r − r k ) = δ k ( r ) is the Dirac delta distribution at the position r k of electron k , and ˆ p k = − i(cid:126) ∇ k is the associated momentum operator. In general, the expectationvalue of any one-electron operator ˆ F can be expressed using Eqs. (3) and (4) as (cid:68) ˆ F (cid:69) ( t ) = (cid:68) Ψ el ( t ) (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) Ψ el ( t ) (cid:69) (12)= (cid:88) λν B † λ ( t ) B ν ( t ) (cid:68) Φ λ (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) Φ ν (cid:69) . (13)Evaluation of the matrix elements (cid:68) Φ λ (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) Φ ν (cid:69) can be done by exploiting the structure of thefunctions { Φ λ , Φ ν } . In the hybrid TDDFT/CIS methodology, these take the form of singly excitedconfigurations, i.e., the truncation of Eq. (4) at the singles level. The matrix elements in the basisof singly excited configurations read (cid:68) Φ λ (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) Φ ν (cid:69) = C ( λ )ref C ( ν )ref (cid:68) φ ref (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) φ ref (cid:69) + (cid:88) bs C ( λ )ref C s ( ν ) b (cid:68) φ ref (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) φ sb (cid:69) + (cid:88) ar C r ( λ ) a C ( ν )ref (cid:68) φ ra (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) φ ref (cid:69) + (cid:88) abrs C r ( λ ) a C s ( ν ) b (cid:68) φ ra (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) φ sb (cid:69) . (14)= C ( λ )ref C ( ν )ref (cid:88) a (cid:68) ϕ a (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) ϕ a (cid:69) + (cid:88) ar C r ( λ ) a C r ( ν ) a (cid:88) a (cid:68) ϕ a (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) ϕ a (cid:69) + (cid:88) bs C ( λ )ref C s ( ν ) b (cid:68) ϕ b (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) ϕ s (cid:69) + (cid:88) ar C r ( λ ) a C ( ν )ref (cid:68) ϕ a (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) ϕ r (cid:69) + (cid:88) ar (cid:54) = s C r ( λ ) a C s ( ν ) a (cid:68) ϕ r (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) ϕ s (cid:69) + (cid:88) a (cid:54) = br C r ( λ ) a C r ( ν ) b (cid:68) ϕ a (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) ϕ b (cid:69) . (15)where a ∈ { , , . . . , a − , r, a + 1 , . . . } denote the occupied spin orbitals of the configuration statefunction | φ ra (cid:105) . Note that we make use of the Slater-Condon rules [8–10] to resolve Eq. (15) in termsof one-electron integrals in the basis of the spin orbitals, ϕ a ( r ). The transition moments betweenspin orbitals are usually computed in the spin-free representation by first integrating over thespin coordinates. Specifically, the expectation value for the electron density requires the followingintegrals (cid:104) ϕ a | ˆ ρ | ϕ b (cid:105) = ρ ab ( r ) = ϕ a ( r ) ϕ b ( r ) . (16)The electronic flux density for a wave packet of the form Eq. (3) can be formulated as j ( r , t ) = 2 i (cid:88) λ<ν Im (cid:104) B † λ ( t ) B ν ( t ) (cid:105) J λν ( r , t ) , (17)which can be calculated by exploiting the CIS structure of the eigenfunctions. The transitionelectronic flux density from state λ to state ν is denoted J λν ( r , t ), which simplifies using the Slater-Condon rules to J λν ( r , t ) = (cid:88) ar (cid:16) C ( λ )ref C r ( ν ) a + C r ( λ ) a C ( ν )ref (cid:17) j ar + (cid:88) a,r (cid:54) = s C r ( λ ) a C s ( ν ) a j rs + (cid:88) a (cid:54) = br C r ( λ ) a C r ( ν ) b j ab (18)6here j ab = − i (cid:16) ϕ a ( r ) (cid:126) ∇ ϕ b ( r ) − ϕ b ( r ) (cid:126) ∇ ϕ a ( r ) (cid:17) (19)are molecular orbital (MO) electronic transition flux densities from MO ϕ a ( r ) to MO ϕ b ( r ).As one of the most widespread bases used in quantum chemistry, we specialize here to spatialMO defined as linear combination of atom-centered orbitals (MO-LCAO for “Molecular Orbital -Linear Combination of Atomic Orbitals”) ϕ a ( r ) = N A (cid:88) A =1 n AO ( A ) (cid:88) i A =1 D ( a ) i A χ i A ( r − R A ) , (20)where D ( a ) i A is the i A th expansion coefficient for MO a . The atomic orbitals χ i A are expressed as afunction of the Cartesian coordinates of one electron r and the spatial coordinates R A of nucleus A . N A labels the number of atoms and n AO ( A ) is the number of atomic orbitals on atom A . Usingthe MO-LCAO ansatz, the transition moments between spin orbitals read (cid:68) ϕ a (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) ϕ b (cid:69) = N A (cid:88) A,B n AO ( A ) (cid:88) i A =1 n AO ( B ) (cid:88) j B =1 D ( a ) i A D ( b ) j B (cid:68) χ i A (cid:12)(cid:12)(cid:12) ˆ F (cid:12)(cid:12)(cid:12) χ j B (cid:69) . (21)The MO-LCAO coefficients D ( a ) iA and the definition of the atomic orbitals can be read directlyfrom the output of standard quantum chemistry program packages. All required derivatives andintegrals in the atomic orbital basis are computed analytically using our Python post-processingtoolbox ORBKIT [11], with which the molecular orbital density (cf. Eq. (16)) and the molec-ular orbital electronic flux density (cf. Eq. (19)) can then be projected on an arbitrary grid.Combining the information in this list with the occupation patterns of the quasi-CI eigenvectorsassociated with the excited states obtained at the LR-TDDFT level of theory, it is possible to cre-ate a library of transition moments between CI-states to be used in the dynamics. Note that thetransition 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 toolsfor the hybrid TDDFT/CIS methodology are implemented, along with various other one-electronquantities, in a recently introduced open-source Python framework detCI@ORBKIT , availableat https://github.com/orbkit/orbkit/ . The program requires a preliminary LR-TDDFT cal-culation using Gaussian-type atom-centered orbitals. There is no restriction for the choice of func-tional. Currently, the code supports the GAMESS [12] and TURBOMOLE [13] formats. Ourprogram then computes matrix elements of one-electron operators, projects them on an arbitrarygrid, and stores them in a library to be used for analyzing the N -electron dynamics. The frame-work detCI@ORBKIT is written in Python, simplifying its portability on different platforms andoffering efficient standard libraries for visualization purposes. Implementation details are given else-where [14]. The dynamics program is not part of the standard implementation and can be performedusing either a user-written code or, e.g., the Matlab WAVEPACKET package [15]. In the present7ork, all dynamical simulations were performed using
GLOCT , an in-house implementation of apropagator for the reduced-density matrix and related quantities [16, 17].
To demonstrate the capabilities of detCI@ORBKIT , we perform the analysis of correlated electrondynamics in two selected molecular systems. First, the charge transfer process in lithium hydrideis studied to benchmark the quality of the TDDFT/CIS description against Full CI results. Theelectron migration in an alizarin dye induced by broadband laser excitation is then used as anexample to demonstrate the scalability of the method.
In lithium hydride, charge migration can be initiated, e.g., by laser excitation from the molecularground state | Ψ g (cid:105) to the first excited state | Ψ e (cid:105) . The charge transfer mechanism can be understoodfrom the analysis of the superposition state | Ψ el ( t ) (cid:105) = 1 √ (cid:16) | Ψ g (cid:105) e − iE g t/ ¯ h + | Ψ e (cid:105) e − i ( E e t/ ¯ h + η ) (cid:17) (22)which leads to the time-dependent one-electron density ρ ( r , t ) = 12 (cid:16) ρ g ( r ) + ρ e ( r ) (cid:17) + ρ ge ( r ) cos(∆ E t/ ¯ h + η ) , (23)where ∆ E = E e − E g , and η is the relative phase. It is set to η = π in the present example.Similarly, the electronic flux density takes the form j ( r , t ) = Im [ J ge ( r )] sin(∆ E t/ ¯ h + η ) , (24)where the transition electronic flux density J ge ( r ) is obtained from Eq. (18). The charge transferdynamics can be thus rationalized in terms of the static transition moment between the two statesinvolved. The electron densities, ρ g ( r ) and ρ e ( r ), of the ground state X Σ + and the charge transferstate A Σ + and the transition density between both, ρ ge ( r ), are computed by combining the MOcontributions obtained from Eq. (16). The Kohn-Sham orbitals, the pseudo-CI eigenvectors, andthe associated LR-TDDFT excitation energies are computed using an aug-cc-pVTZ at the B2-PLYPlevel of theory, as implemented in TURBOMOLE. [13] The character of the charge transfer state isfound to be dominated by the HOMO-LUMO transition (see Fig. 1 (right side)). This is in goodagreement with the character determined from CIS and Full CI calculations, both performed withPSI4 [18] using the identical basis set. The corresponding frontier orbitals from the Hartree-Fockreference 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 thisdifference has only a marginal influence on the electron dynamics.8 great advantage of LR-TDDFT over CIS is the improved energetic description of the excitedstates at virtually the same computational cost. This is demonstrated by the good agreement ofthe excitation energies at the B2-PLYP level of theory (∆ E = 3 .
51 eV) with the Full CI reference(∆ E = 3 .
56 eV), as compared to the rather poor value for the truncated CIS expansion (∆ E =4 .
04 eV). Since the excitation is dominated by single excitations in all cases, the transition energieswill affect mostly the timescale of the dynamics. Further, considering the similarities between theMOs involved in the dominant transition, a similar dynamical behavior is to be expected. This isindeed what is observed in Fig. 2, where the left panels report the flux density j ( r , t ) at time t = τ / y ( r , t ) computed fromthe one-electron density, Eq. (23). In the top right panel, the benchmark Full CI calculation revealsthat 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 recognizedaround 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 aroundthe hydrogen for the two single electron approaches.The electronic flux density maps are depicted as streamlines in the left panels of Fig. 2, wherethe blue shades indicate its magnitude. The charge is seen to flow from the polarizable hydrogenaround and towards the back of the hard lithium. The same large vortex around the lithiumion is observed for the Full CI benchmark, the CIS, and the hybrid TDDFT/CIS approach. Themain quantitative differences between the methods are located in the low-density regions, e.g., at x > a , where the Full CI benchmark predicts a flux almost parallel to the molecular axis. Thecritical point (at x ∼ a ) between the lithium and hydrogen atoms also appears to be slightlyshifted to the right at the CIS and TDDFT/CIS level of theory. In general, it can be said that bothsingle electron excitation ansatz yield a very similar picture of the electron difference density andthe associated flux density. This conclusion is likely to hold for all dynamical processes involving N -electron eigenstates dominated by a single excitation character. In this second example, we demonstrate the computational scalability of the hybrid TDDFT/CISapproach to analyze the correlated electron dynamics for more extended molecular systems. Thenecessity 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 thesingles excitation represents a simple, intuitive, and computationally cheap approach to computequalitatively correct excited electronic states. [19,20] However, it yields inaccurate vertical excitationenergies from the ground state. [21] As advocated in the theory section above, a suitable alternativeis LR-TDDFT [22], which usually provides better energetic description than CIS while retaining thesame 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 domi-9ated by a single excitation character. In addition, it benefits from the versatility and continuousimprovement of density functionals. Proper treatment of the excited states strongly depends on theappropriate 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–29] Fortunately,extensive experience has been accumulated over the years concerning the applicatibility of eachfunctional in specific chemical contexts. For example, it is known that non-local exchange improvesgreatly the description of charge transfer states [20].To show the scalability of the hybrid TDDFT/CIS formalism for the analysis of correlatedelectron dynamics, we initiate a photoinduced ultrafast charge migration process in alizarin. Thisorganic chromophore is used as a π -conjugated photosensitizer in dye-sensitized solar cells. Prior tothe dynamical simulation, the electronic and optical properties of alizarin are determined by meansof LR-TDDFT. Therefore, we perform a TURBOMOLE [13] calculation with the B3LYP [30] hybridfunctional and the def2-SVP basis set [31,32] at the equilibrium geometry of alizarin. This setup hasbeen previously proven to yield accurate results for the electronic and spectroscopic properties ofsuch systems. [4, 33–35] The computed absorption spectrum is depicted in Fig. 3(a) along with theexperimentally observed absorption band of free alizarin (dashed black line). The good agreementfor the first absorption band between theory (437 nm) and experiment (431 nm) underlines thesuitability of TDDFT to model the electronic spectra of medium-sized organic molecules dominatedby single excitation character. It is important to recognize that the second absorption band iscomposed 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 broad-band 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 sin -shaped pulses with a duration 19 fs is constructed. The pulse is adjusted to the parametersof a realistic experimental laser field used in similar investigations to initiate, e.g., ultrafast photoin-duced processes in alizarin-TiO solar cells. [36, 37] The resulting electric field is shown in the insetof Fig. 3(b). The laser excitation is followed by a 20 fs field-free propagation. The time-evolution ofthe N -electron wave packet (cf. Eqs. (1) and (3)) is accomplished using an adaptive Runge-Kuttaalgorithm in the interaction picture. The methodology and implementation details are describedelsewhere. [16, 38, 39] Fig. 3(b) shows the evolution of the state populations and the applied laserfield in the inset. As it is often the case for molecules in strong fields, the population dynamics isvery intricate while the laser is on, in part due to important polarization effects and in part due tothe number of states that are excited by the broadband laser. To account for the electronic responseof the system to the laser field, 25 eigenstates are incorporated in the simulation. After the laserexcitation, only twelve states are significantly populated ( P λ > . N -electron wave packet in the pseudo-CI eigenvector basis. This10an become computationally tedious, since the number of Slater determinants in the wave functionexpansion increases with the number of occupied and unoccupied orbitals. For alizarin in thecurrent basis set, the 62 occupied MOs times 248 virtual MOs correspond to 15 376 determinantsfor each of the excited states. Recalling that reconstruction of the flux density, Eq. (18), requirescombining one-electron integrals of all orbitals pairwise, for each pair of eigenstates, this amountsto a tremendous computational task. However, inclusion of the complete set of determinants isnot 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 implementationand can be controlled by a user-defined convergence threshold.In order to examine the influence of truncating the complete set of determinants to a physicallymeaningful subset in the wavefunction ansatz, we choose to use a threshold for the expansioncoefficients of the quasi-CI wave functions, Eq. (4). All determinants with a contribution under thisvalue are neglected in the evaluation of the one-electron density and the associated flux density.We define three different thresholds, (cid:12)(cid:12)(cid:12) C r, ( λ ) a (cid:12)(cid:12)(cid:12) > { − , − , − } . To illustrate the computationalsavings that can be expected, the thirteenth excited state is selected as an example, since it is themost populated state after the broadband laser excitation. The three chosen threshold values retain926, 31, and 7 Slater determinants in the wave function expansion, respectively. For consistencycheck, the pruned expansions recover 99.94 %, 99.31 %, and 95.92 % of the norm of the thirteenthexcited state, respectively. These numbers are similar for other excited states. After excludingthe contributions lying under the respective threshold, the remaining coefficients are renormalizedusing a modified Gram-Schmidt algorithm. Regarding the computational effort, the reduction ofdeterminants means a drastic decrease of computational steps. For example, the calculation ofthe transition electronic flux density j ( r ) between the fifth and thirteenth excited state requires N × N = 573 ×
926 = 530 598 determinant combinations for the more strict threshold of (cid:12)(cid:12)(cid:12) C r,λa (cid:12)(cid:12)(cid:12) > − . This number reduces to N × N = 3 × (cid:12)(cid:12)(cid:12) C r,λa (cid:12)(cid:12)(cid:12) > − . This simple constraint thus confers great scalability to the methodpresented here.To assess the quality of this approximation, the electronic flux densities j ( r , t ) reconstructedusing the different thresholds are illustrated in Fig. 4(a),(c),(d) at a characteristic point in timeafter the laser-pulse excitation (here, t = 25 . ∂ρ ( r , t ) /∂t , is displayed additionally for the tight threshold, (cid:12)(cid:12)(cid:12) C r, ( λ ) a (cid:12)(cid:12)(cid:12) > − , tofacilitate the interpretation of the flux densities. For all three wave packet expansions, j ( r , t ) showsnearly identical qualitative features. These include: (1) an electron flow along the bonds mediated11y the π -system, (2) an anti-clockwise flux at the outer right ring, and (3) a charge migration fromthe outer right ring and the top part of the left ring to the lower hydroxyl group. Due to the cyclicnature 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 alllevels of analysis, the very important numerical savings at the crudest level of approximation confergreat 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 thelower hydroxyl group to the neighboring carboxyl group is not fully reproduced with the two largerthresholds (cf. Fig. 4(c) and (d)). This corresponds to a through-space charge transfer and ismediated by a large number of small contributions in the pseudo-CI eigenstate basis. While themajor phenomenological characteristics of the correlated electron dynamics are still captured, someminor mechanistic information is lost by reducing the wave functions to their dominant determinantcontributions.To understand the electron dynamics, we extend its analysis to the time-dependent electronicyield. It is defined as the difference between the electron density at a given time and the electrondensity at t = 0 fs, integrated over a given volume, y A ( t ) = (cid:90) A d r (cid:16) ρ ( r , t ) − ρ ( r , (cid:17) (25)In Fig. 5, the difference densities projected on right and left aromatic ring of alizarin are reportedfor times after the laser has been switched off. These reveal intricate synchronous fluctuations ofelectron density between both rings. The slight asymmetry of the electron redistribution correlateswith the asymmetric substituents on the two outer rings: the fluctuation amplitudes in the left ringare larger than in the unsubstituted, rightmost aromatic ring. In Fig. 6, the electronic flux densityis plotted for selected characteristic points in order to unveil the mechanism of the charge migrationbetween both rings. The times associated with these snapshots are marked as vertical gray lines inFig. 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 magneticfield. In the next snapshot (panel (b), t =25.3 fs), the yields in the left and right aromatic ringsare 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 thebonds of the right ring, passing by the bottom carbonyl group to bottom hydroxyl group belowthe left aromatic ring, and back to the central carbonyl unit via a through space mechanism. Thisis a strong indication of the presence of a hydrogen bond, which will lead to a rapid tautomerichydrogen transfer. At the third time step ( t =25.9 fs), the electronic yield is maximal in the rightring and the electron flux density rotates clockwise. Contrary to panel (a), the electron flow is notas strongly localized inside the ring, probably due to the asymmetry of the ring substituents. In thelast snapshot ( t =26.4 fs), the density migrates back from the leftmost to the rightmost aromaticring along a different path: the electrons retain a clockwise orientation in the left ring, mostly12owing from the bottom hydroxyl to the top central carbonyl group. A marginal amount flowsfrom the bottom left hydroxy group over the hydrogen bond to the bottom carbonyl group andinto rightmost aromatic ring. For the complete picture of the dynamics after the laser excitation, ashort film of the charge migration process is made available online in the Supporting Information. In this paper, we have introduced a novel procedure to analyze and visualize many-electron dynamicsfrom 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 basisof 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 retainsthe simple character of coupled electron-hole pairs, which facilitates its interpretation in terms ofconfiguration states while keeping the size of the basis relatively small. This renders the hybridmethod amenable to large systems – in fact, any system that can be tackled and described accuratelyusing linear-response TDDFT.From the evolution of a TDDFT/CIS wave packet, we showed how it is possible to rationalizethe N -electron dynamics of a system in terms of transition moments of various one-electron op-erators. These include difference electronic densities, transient electronic flux density maps, andtransition dipole moments, which we have implemented in a Python toolbox for postprocessingmulti-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 projectedon an arbitrary grid. The required information (molecular orbitals, structure, Gaussian atomic ba-sis, pseudo-CI coefficients) can be directly extracted from the output of various quantum chemistryprogram packages. These are then used to reconstruct quantities that help understanding the flowof electrons in molecules out of equilibrium.We first applied the hybrid TDDFT/CIS scheme to a test system, the charge transfer in lithiumhydride, in order to benchmark the quality of the ansatz against standard CIS and Full CI referencesimulations. The results demonstrated that a good choice of the functional improves mostly theenergetic description of the charge transfer state, which can be brought close to the Full CI bench-mark. On the other hand, the pseudo-CI basis retains a single excitation character, which is foundto be similar to the standard CIS reference. Both ansatz agree semi-quantitatively with the higherlevel 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 chromophorefor dye-sensitized solar cells demonstrated the scalability of the method and the versatility of thenew toolkit. In particular, it was found that the main features of the electron flow mechanismcan be recovered using a stringent basis pruning strategy, in which each pseudo-CI eigenstate isrepresented using only a few dominant configurations. By doing so, marginal features of the many-13lectron dynamics involving, e.g., electron flow through hydrogen bonds, may be lost. Because ofthe favorable scaling of the method even at tighter convergence Thresholds, it is expected to beapplicable to a large number of medium-sized molecules. This will help to understand electronmigration processes in various fields, such as nano electronics and light-harvesting applications.
The authors gratefully thank the Scientific Computing Services Unit of the Zentraleinrichtung f¨urDatenverarbeitung at Freie Universt¨at Berlin for allocation of computer time, and Hans-ChristianHege for providing the ZIBAmira visualization program. The funding of the Deutsche Forschungs-gemeinschaft (project TR1109/2-1) and of the Elsa-Neumann foundation of the Land Berlin arealso acknowledged.
Correlated Electron Dynamics, Time-dependent Density Functional Theory, Electronic Flux Den-sity, Electron Density, Electronic Current Density14 eferences [1] M. Dantus, M. J. Rosker, and A. H. Zewail. Real-time femtosecond probing of transition statesin chemical reactions.
J. Chem. Phys. , 87:2395, 1987.[2] T. S. Rose, M. J. Rosker, and A. H. Zewail. Femtosecond real-time observation of wave packetoscillations (resonance) in dissociation reactions.
J. Chem. Phys. , 88:6672, 1988.[3] J. A. Sonk, M. Caricato, and H. B. Schlegel. Td-ci simulation of the electronic optical responseof molecules in intense fields: Comparison of RPA, CIS, CIS(D), and EOM-CCSD.
J. Phys.Chem. A , 115:4678–4690, 2011.[4] G. Hermann and J. C. Tremblay. Ultrafast photoelectron migration in dye-sensitized solarcells: influence of the binding mode and many-body interactions.
J. Chem. Phys. , 145:174704,2016.[5] S. Klinkusch and J. C. Tremblay. Resolution-of-identity stochastic time-dependent con-figuration interaction for dissipative electron dynamics in strong fields.
J. Chem. Phys. ,144(18):184108, 2016.[6] S. Pirhadi, J. Sunseri, and D. R. Koes. Open source molecular modeling.
J. Mol. Graph.Model. , 69:127, 2016.[7] E. Schr¨odinger. Quantisierung als Eigenwertproblem.
Ann. Phys. (Leipzig) , 81:109, 1926.[8] J. C. Slater. The theory of complex spectra.
Phys. Rev. , 34(10):1293, 1929.[9] E. U. Condon. The theory of complex spectra.
Phys. Rev. , 36(7):1121, 1930.[10] J. C. Slater. Molecular energy levels and valence bonds.
Phys. Rev. , 38:1109, Sep 1931.[11] G. Hermann, V. Pohl, J. C. Tremblay, B. Paulus, H.-C. Hege, and A. Schild.
ORBKIT : Amodular python toolbox for cross-platform postprocessing of quantum chemical wavefunctiondata.
J. Comput. Chem. , 37(16):1511, 2016.[12] Michael W. Schmidt, Kim K. Baldridge, Jerry A. Boatz, Steven T. Elbert, Mark S. Gordon,Jan H. Jensen, Shiro Koseki, Nikita Matsunaga, Kiet A. Nguyen, Shujun Su, Theresa L.Windus, Michel Dupuis, and John A. Montgomery. General atomic and molecular electronicstructure system.
J. Comput. Chem. N -Electron Dynamics: I. Multi-Determinantal Wave Functions. ArXiv e-prints arXiv:1701.06885[physics.chem-ph] , January 2017.[15] Ulf Lorenz B. Schmidt. WavePacket 5.2, : A matlab program package for quantum-mechanicalwavepacket and density propagation and time-dependent spectroscopy., 2016. Available viahttp://sourceforge.net/p/wavepacket/matlab.[16] J. C. Tremblay, T. Klamroth, and P. Saalfrank. Time-dependent configuration-interactioncalculations of laser-driven dynamics in presence of dissipation.
J. Chem. Phys. , 129(8):084302,2008.[17] J. C. Tremblay and P. Saalfrank. Guided locally-optimal control of quantum dynamics indissipative environments.
Phys. Rev. A , 78:063408:1–9, 2008.[18] J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evangelista, J. T.Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, N. J. Russ, M. L. Leininger,C. L. Janssen, E. T. Seidl, W. D. Allen, H. F. Schaefer, R. A. King, E. F. Valeev, C. D.Sherrill, and T. D. Crawford. Psi4: an open-source ab initio electronic structure program.
WIREs Comput. Mol. Sci. , 2(4):556, 2012.[19] J. B. Foresman, M. Head-Gordon, J. A. Pople, and M. J. Frisch. Toward a systematic molecularorbital theory for excited states.
J. Phys. Chem. , 96(1):135, 1992.[20] A. Dreuw, J. L. Weisman, and M. Head-Gordon. Long-range charge-transfer excited statesin time-dependent density functional theory require non-local exchange.
J. Chem. Phys. ,119(6):2943, 2003.[21] A. Dreuw and M. Head-Gordon. Single-reference ab initio methods for the calculation of excitedstates of large molecules.
Chem. Rev. , 105(11):4009, 2005.[22] E. K. U. Gross and W. Kohn. Time-dependent density-functional theory.
Adv. QuantumChem. , 21:255, 1990.[23] M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub. Molecular excitation energies tohigh-lying bound states from time-dependent density-functional response theory: Characteri-zation and correction of the time-dependent local density approximation ionization threshold.
J. Chem. Phys. , 108(11):4439, 1998.[24] H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao. A long-range correction scheme for generalized-gradient-approximation exchange functionals.
J. Chem. Phys. , 115(8):3540, 2001.[25] J. Heyd, G. E. Scuseria, and M. Ernzerhof. Hybrid functionals based on a screened coulombpotential.
J. Chem. Phys. , 118(18):8207, 2003.1626] Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai, and K. Hirao. A long-range-corrected time-dependent density functional theory.
J. Chem. Phys. , 120(18):8425, 2004.[27] M. A. Rohrdanz and J. M. Herbert. Simultaneous benchmarking of ground-and excited-stateproperties with long-range-corrected density functional theory.
J. Chem. Phys. , 129(3):034107,2008.[28] M. A. Rohrdanz, K. M. Martins, and J. M. Herbert. A long-range-corrected density functionalthat performs well for both ground-state properties and time-dependent density functional the-ory excitation energies, including charge-transfer excited states.
J. Chem. Phys. , 130(5):054112,2009.[29] D. Jacquemin, V. Wathelet, E. A. Perpete, and C. Adamo. Extensive TD-DFT benchmark:singlet-excited states of organic molecules.
J. Chem. Theory Comput. , 5(9):2420, 2009.[30] A. D. Becke. Density-functional exchange-energy approximation with correct asymptotic be-havior.
Phys. Rev. A , 38(6):3098, 1988.[31] A. Sch¨afer, H. Horn, and R. Ahlrichs. Fully optimized contracted gaussian basis sets for atomsLi to Kr.
J. Chem. Phys. , 97(4):2571, 1992.[32] F. Weigend and R. Ahlrichs. Balanced basis sets of split valence, triple zeta valence andquadruple zeta valence quality for H to Rn: design and assessment of accuracy.
Phys. Chem.Chem. Phys. , 7(18):3297, 2005.[33] W. R. Duncan and O. V. Prezhdo. Electronic structure and spectra of catechol and alizarin inthe gas phase and attached to titanium.
J. Phys. Chem. B , 109(1):365, 2005.[34] R. S´anchez-de Armas, J. Oviedo L´opez, M. A. San-Miguel, J. F. Sanz, P. Ordej´on, andM. Pruneda. Real-time TD-DFT simulations in dye sensitized solar cells: the electronic ab-sorption spectrum of alizarin supported on TiO nanoclusters. J. Chem. Theory Comput. ,6(9):2856, 2010.[35] T. Gomez, G. Hermann, X. Zarate, J. P´erez-Torres, and J. C. Tremblay. Imaging the ultrafastphotoelectron transfer process in alizarin-TiO . Molecules , 20(8):13830, Jul 2015.[36] R. Huber, J.-E. Moser, M. Gr¨atzel, and J. Wachtveitl. Real-time observation of photoinducedadiabatic electron transfer in strongly coupled dye/semiconductor colloidal systems with a 6 fstime constant.
J. Phys. Chem. B , 106(25):6494, 2002.[37] L. Dworak, V. V. Matylitsky, and J. Wachtveitl. Ultrafast photoinduced processes in alizarin-sensitized metal oxide mesoporous films.
ChemPhysChem , 10(2):384, 2009.1738] J. C. Tremblay and T. Carrington Jr. Using preconditioned adaptive step size Runge-Kuttamethods for solving the time-dependent Schr¨odinger equation.
J. Chem. Phys. , 121(23):11535,2004.[39] J. C. Tremblay, S. Klinkusch, T. Klamroth, and P. Saalfrank. Dissipative many-electron dy-namics of ionizing systems.
J. Chem. Phys. , 134(4):044311, 2011.18igure 1: Two-dimensional projections of the dominant orbitals involved in the excitation of thefirst charge transfer state in LiH. The HOMO is found to be similar at Hartree-Fock (left) and atthe B2-PLYP (right) levels of theory. The LUMO obtained from B2-PLYP is more delocalized andfound at significantly higher energy than the comprable Hartree-Fock orbital.19igure 2: Comparison of the electron transfer process in LiH at time t = τ /
4, obtained usingdifferent electronic structure methods (a) Full CI; (b) CI Singles; (c) hybrid TDDFT/CIS usingthe B2-PLYP functional. Left panels: the black arrows show the electronic flux density, j ( r , t ),represented as streamlines. The color bar indicates the regions of high flux density (dark blue) andlow flux density (white) in units of E h / ¯ ha . Right panels: the difference density in units of a − ,Eq. (9), illustrates that an electron is transferred from the hydrogen (density depletion in blue)to the lithium ion (density increase in red). All qualitative features are in good agreement for allmethods. 20a)
200 300 400 500 hc/ ∆ E [nm]0 . . . . . . . O s c ill a t o r S tr e n g t h (b)Figure 3: (a) Simulated optical spectrum of alizarin obtained from linear-response TDDFT usingthe B3LYP functional and a def2-SVP basis set. The vertical gray lines represent the oscillatorstrengths of specific transitions. The broadened spectrum (solid blue line) is constructed usingGaussian functions with width σ = 25 nm. The position of the experimental band of the freealizarin at 431 nm is marked by a dashed black line. (b) Population evolution of the ultrafast N -electron dynamics in alizarin driven by broadband laser excitation. The inset shows the timeevolution of the laser pulse. All states between 200 nm and 500 nm are excited from the groundstate. 21a) (b)(c) (d)Figure 4: Comparison of the electronic flux density j ( r , t ) calculated on the basis of a pseudo-CIeigenfunctions with three different cutoff thresholds, i.e., (a) | C r, ( λ ) a | > − , (c) | C r, ( λ ) a | > − ,and (d) | C r, ( λ ) a | > − . The vector arrows are colored according to their magnitude. In addition,the electron flow ∂ρ ( r , t ) /∂t calculated with the stringent cutoff threshold is depicted as a contourplot. The negative and positive isosurfaces are colored gray and blue with a isosurface value of ± · − a − . The time ( t = 25 . t [fs] − . − . − . . . y A t t t t Figure 5: Time-evolution of the yield y A projected on the left (black) and right (blue) aromaticring of the alizarin dye after the broadband laser excitation. The vertical gray lines representcharacteristic time points for which the time-dependent electronic flux densities j ( r , t ) are evaluatedin Fig. 6. 23a) (b)(c) (d)Figure 6: Vector plots of the electronic flux density j ( r , t, t