A quantum dynamics method for excited electrons in molecular aggregate system using a group diabatic Fock matrix
Abstract
We introduce a practical calculation scheme for the description of excited electron dynamics in molecular aggregate systems within a local group diabatic Fock representation. This scheme makes it easy to analyze the interacting timedependent excitation of local sites in complex systems. In addition, lightelectron couplings are considered. The present scheme is intended for investigations on the migration dynamics of excited electrons in lightinduced energy transfer systems. The scheme was applied to two systems: a naphthalene(NPTL)tetracyanoethylene(TCNE) dimer and a 20mer circle of ethylene molecules. Through local group analyses of the dynamical electrons, we obtained an intuitive understanding of the electron transfers between the monomers.
I Introduction
The interaction of light with molecular systems is a fundamental factor for energy source in nature. The photo chemistry as a representative process involved with this type of interaction is characterized by excited electrons in molecular systems. Especially, the dynamics of excited electrons in aggregated molecular systems play a significant role in energy transfer phenomena. KuhnCET (); Scholesrevexcitons (); SchifferRevPCET (); SpanoEET () Electron, radical and exciton transfers are important topic in a research field treating excited electron dynamics in molecular aggregates. The electron transfer is characterized by inter or intra molecular charge migration. The radical electron is well described by unpaired singly occupied electrons in a natural orbital representation. The exciton transfer is caused by interactions of excitons created in coherently coupled monomer group as a holeelectron pair often accompanied with their delocalization. From the viewpoint of real time dynamics in a light energy conversion, several pioneering investigations have been reported in the fields of photo synthesis, solar cells, and photo chemical reactions. artps (); rtexppcbm (); rtexpGraetzel (); envassistqweet ()
Theories for the treatment of realtime electron dynamics in complex excited states accompanied by the changes in chemical bonding, have been developed and applied, for example, realtime time dependent density functional theory, time dependent configuration interaction, and configuration state function schemes. RGTDDFT (); octpusrttddft (); gceed (); PCCPeldtKTTY (); NonadiabaticityChemRev (); ChemicaltheorybeyondBOparadigm () These theories provide an excellent description of electron dynamics in total systems. However, there is no explicit scheme for the analysis of the timedependent interaction dynamics of excited electrons in a complex system with using a reliable level of ab initio electronic structure calculation. The schemes described are successful but timeconsuming, so a general theoretical tool allowing the compact description and investigation of the electron dynamics in any aggregated molecular would be useful.
In the present work, we propose a practical theory for describing the excited electron dynamics in aggregated molecular systems driven by the interactions between subsystems and by lightmatter interactions.
This treatment enables us to treat arbitrary types of initial local excitation in a complex system. We have used the wellknown group diabatic matrix representation of Fock operator. Lan (); Gao (); Shi (); ThossTiO2 () We call this the GDF model. A skeleton Fock matrix of the total system is represented by a set of localized orbitals of the sub groups after applying a Löwdin orthogonalizaion of the atomic orbitals. Szabo () In this article, we show that this provides a useful scheme for studying the electron dynamics accompanied by internally/externally induced excitation. Combined with lightmatter couplings and the skeleton Fock matrix mentioned above, we can extract the essential features of the electron migration dynamics at an early stage after any type of excitation. Thus, our scheme can, potentially, provide meaningful information for the study of the transfer of lightenergy to chemical functionality via molecular aggregate systems. We show the efficiency of the introduced theoretical method by two example applications, a typical donoracceptor dimer system and an aggregate of circularly arranged monomers under external continuum radiation fields. The features found in the dynamics of the excited electrons in these aggregated systems validate the efficiency of the present scheme.
The contents of this report are as follows: In Sect. II, we review current theoretical methods concerning electron dynamics using the equations of motions and explain the need to introduce a group interaction representation to the electron dynamics theory. Then we explain the details of the construction of the GDF representation. Further details of the equation of motion of the density matrix in this representation are also given in this section. In Sect. III, we present numerical applications of the present method by using two types of aggregated systems. In the last section IV, we provide concluding remarks.
Ii Theoretical method
ii.1 Computational theories for the study electron dynamics
Here, we review the theories of electron dynamics in molecular systems and we explain the need for our newly devised theoretical method in this article. Although the interplay between the nuclear and electronic wave packets is an important topic, PSANBlaser (); encoinsJPCA (); BohmianSuzuki2016 (); FactorizationENGross () we have not discussed this to allow us to focus on our electron dynamics description scheme.
There are several types of practical theories for the description of electron dynamics involving excited states in molecular systems. These theories are based on the electronic structure calculations. For brevity, they can be roughly categorized into two schemes: those directly utilizing (a) a time dependent multielectron wavefunction as a timedependent linear combination of configuration state function constructed from a molecular orbital set and (b) a timedependent electron density matrix made from independent orbitals propagating at the electron timescale.
First, in the following two subsections, we review separately these two schemes, and in the last subsection, we address the differences with a related scheme for the treatment of the interaction dynamics in excited aggregated systems.
ii.1.1 (a) Schemes based on multistate multielectron wavefunctions
The timedependent configuration interaction theory (TDCI) including the time dependent configuration state function (TDCSF) SaalfrankJCP2007TDCI (); TDCSFAmano (); NonadiabaticityChemRev (); ChemicaltheorybeyondBOparadigm () and multi configurational timedependent HartreeFock (MCTDHF) schemes mctdhfscrinzi (); mctdhfkatokono () are representatives of the class (a).
In the TDCI scheme, the propagation of the electron wave packet obeys a snapshot equation of motion at each molecular configuration. Here, an electron wave packet is defined as the superposition of multielectronic adiabatic state functions of a molecular system. Only the configuration interaction (CI) coefficients are propagated in the time domain while the electron orbitals for constructing the configuration basis functions are determined in a static manner at each molecular geometry. In cases considering the molecular structure dynamics, the orbitals uniformity as well as resultant uniformity of CSF made from orbitals with changes in the molecular geometry, are needed for consistency between the timedependent CSFs and their coefficients with dynamics. This includes the solution responsible for the orbital character exchange with molecular motion. The type depends on the variation of the methods for the compact, smooth, and correct description of the electron dynamics. TDCSFAmano (); NonadiabaticityChemRev (); TDCSFKunisadaUshiyama ()
Unlike TDCI, where only the coefficients of the CI functions or CSFs propagate in the time domain, in the MCTDHF scheme, which is a Fermion system counterpart of the multi configuration time dependent Hartree (MCTDH) method, bookMeyerMCTDH () both the configuration coefficients and orbital parameters are simultaneously propagated according to the equation of motion derived from the application of the timedependent variational principle to the quantum action. Here, timedependent variational principle means the equation of motions for the parameters determining a time dependent quantum state is derived by the stationalization of the quantum action as a functional of time dependent quantum state, , under the constraint of normalization condition with respect to the variation of where denotes Hamiltonian operator of a system under consideration. Because of this substantially instantaneous optimization of the time dependent orbitals by using the effective oneelectron selfconsistent field, in MCTDHFtype of schemes, we can save the dimensions of the multielectron configuration functions. This approach is usually applied to highaccuracy studies of field induced electron dynamics in small molecular systems with high accuracy. As practical schemes accompanied by the modelling of a timedependent active orbital space in a sophisticated manner, the timedependent complete active space selfconsistent field (TDCASSCF) tdcascfsato () theory and the restricted active space counterpart (TDRASSCF) tdrasscfMadsen () were developed in a similar manner to that of MCTDHF theory.
One advantage of the schemes belonging to (a) is that they offer a picture of the electron wave packet dynamics in terms of multiple multielectron adiabatic or pseudodiabatic states. For example, this aspect allows the theoretical analysis of nonadiabatic and optical transitions between multi electronic adiabatic states. The methods categorized into (a) can principally reap the benefit of high level electronic correlation methods including postHartreeFock theories. TDCASCISDNest (); TDDMRG ()
Applications with using (a) include a fieldinduced electron dynamics at the attosecond timescale, Diboranelaser (); PCCPeldtKTTY () nonadiabatic electron wave packets in a protoncoupled electron transfer processes, UshiyamaACIE2006PCET (); nagashima2012 (); TDCSFpcetyamataka () the characterization of electron dynamics in highly quasidegenerate excited states b12nad (), and ultrafast photoionization process. MuCurdyMCTDHFapplionization (); TDCSFionizationmatsutaka () Despite the successes in analyzing electron dynamics at the attosecond scale, the system size is inevitably limited by the computational cost associated with the complex mathematical ingredients of these methods, thus reducing the ease of their use.
ii.1.2 (b) Schemes based on the electron density matrix
The class (b) includes the time dependent HartreeFock (TDHF) and, KulanderPRArtTDHF () realtime timedependent density functional theories (RTTDDFT). RGTDDFT (); octpusrttddft (); gceed (); MengKaxirasJCP2004rttddftAO () Here, the electron density operators and associated matrices play central roles in the description of realtime electron dynamics. In particular, compared to the class (a) RTTDDFT provides a way to treat large molecular systems by considering the electronic correlation at a moderate level using effective approximations, resulting in a reasonable computational cost.
There are two types for RTTDDFT. One utilizes timedependent mutually independent orthogonal orbitals obeying the timedependent KohnSham (TDKS) equation, RubioJCP2004propagationtdks () where a timedependent one electron field is determined by the electron density matrix constructed via the ground state occupation of the TDKS orbitals in the same way as that of density functional theory (DFT). Note that the determination of the occupied orbitals is equivalent to that of the oneelectron density matrix. In this orbital scheme, all the occupied orbitals are directly propagated at the electron timescale. The equation of motion of the th TDKS orbital, , is written as
(1) 
where is an orbital label and runs over all occupied orbitals. The KohnSham Fock operator, , depends on the occupied independent KohnSham orbitals. This is expressed by Eq. (1) as a set of functions in the parentheses , attached to . Here means whole the set of occupied TDKS orbitals under consideration.
The other scheme of RTTDDFT is based on the direct time propagation of the electron density matrix without any propagation of the KS orbitals. VoorithrrtTDDFTdensitymat () The KS orbitals are determined through a selfconsistent field calculation that is carried out once at the initial simulation time and constructs the initial electron density matrix needed for a time propagation. Formally, the time propagation of density matrix, , within any basis set representation obeys the Liouville von Neumann equation, including the one electron effective Hamiltonian, which is constructed as the sum of the timedependent Fock matrix, , and lightelectron coupling Hamiltonian, ,
(2) 
Here, . Note that the Fock matrix depends on the timedependent density matrix at each electron propagation timestep. Therefore, the updating the constructed Fock matrix is the dominant part of the computation. Usually some approximations with respect to the update of the Fock matrix and propagation schemes in this nonlinear type equation of motion are employed to reduce computational cost.
These two schemes using orbitals and densities are equivalent to each other within the TDDFT framework under the adiabatic approximation for electron functionals named as timeinlocal one and orbital occupation manner for the time dependent KS orbitals as referential independent orbitals. bookfundamentalTDDFTGross () In the main part of this article, we use the density propagation scheme as explained in this subsection but with an extra approximation.
ii.1.3 Group analysis of electron dynamics in an aggregated system
In previous studies, a time dependent interaction analysis was not sufficient to investigate the excited electron migration dynamics in molecular aggregate systems. Below, we describe two theoretical works related to our approach and our aims.
The first approach uses the molecular orbital (MO)CI Quantum Master Equation (MOQME) method MOQME () developed by Kishi and Nakano. This theory can treat the electron dynamics in an excited molecular aggregate by using a quantum master equation and density matrix within a CI representation. This has been successfully applied to donoracceptor systems and has revealed the mechanism of exciton recurrence motions that are dependent on an external laser fields. The differences between this method and ours are threefold. The first difference concerns the type of ab initio calculation. Kishi and Nakano used density matrices associated with the CIS multi electron wavefunctions. In contrast, our scheme is a Kohn Sham Fock based method that includes effective dynamic correlation. The second difference is that Kishi and Nakano used the BornMarkov approximation based on the weak interaction of monomers in the aggregated system; we did not use this type of approximation to avoid the inherent loss of information, which would limit the future extension of our method. The third difference is that our scheme can use any group separation for group local diabatization in the representation of electron property matrices, including the density matrix, which allows the explicit analysis of group interactions in a timedependent manner.
The second approaches are the study of electron transfer dynamics in a test molecular donoracceptor pair for use as a solar material and a realistic model of photo energy conversion material by using the GDF Hamiltonian and the advanced version of MCTDH method (multilayer MCTDH) by Thoss ThossTiO2 () and Xie et al. Lan (), respectively. The former author, Thoss, is one of the founders of the multilayer MCTDH MLMCTDH () scheme which made breakthrough in a higher dimensional quantum dynamics study and also made intensive researches by using this block diagonal approach.
Thoss et al. examined the photoinduced electrontransfer process in the alizarinTiO2 system as a dyesemiconductor system by using the vibronic model Hamiltonian obtained by the first principle calculation and the block diagonalization approach. They found an interesting feature that an electron injection process proceed in a femtosecond time scale, which is assisted by a significant electronic coherence.
In the work by Xie et al., they constructed a fourstate electronic Hamiltonian using a GDF representation and DFT calculations and then investigated the multielectronicstate nuclear wave packet dynamics using 4 electronic states and 246 vibrational modes. They carefully chose the electronic states responsible for the charge transfer dynamics and, by increasing the number of vibrational modes, obtained valuable information concerning the effect of the vibrational motion on the electron dynamics. Although their work is excellent and partly analyzed the electronic properties as diabatic state populations by using the nuclear wave packet dynamics, it remains challenging to examine the associated electronic properties based on a localized representation, as well as, for example, unpaired electrons and bond order properties, among others. For example, in the cases with external laser fields, many more electronic states should be considered. However, this is difficult to treat using the approaches mentioned above.
A transparent as well as convenient group analysis scheme for describing electron dynamics is required, not only for the investigation of the unresolved excited electron affinities between or within constituent molecules. Scholesrevexcitons ()
In the following parts of this article, a description of our theoretical methodology for this purpose is given, and its validity is verified through application to two types of aggregated systems.
ii.2 Group diabatic representation
On constructing a group diabatic representation of the effective electronic Hamiltonian for the examination of the excited electron dynamics in a general molecular aggregated system, we employed the procedures applied in the articles of Lan, Gao, Shi and Thoss. Lan (); Gao (); Shi (); ThossTiO2 () However, here, we have added a density matrix time propagation scheme and lightelectron couplings to their formulation. Below we summarize this briefly.
ii.2.1 Fock matrix with orthogonalized Löwdin atomic orbitals
We employ the representation matrix of the Fock operator in terms of the wellknown Löwdin orthogonalized atomic basis function, Szabo ()
(3) 
where the orthogonalized Löwdin atomic orbitals(AO) are expressed by
(4) 
with being the AO overlap matrix element and being the original AOs. The original AO representation matrix of the Fock operator, , is responsible for the electron density of arbitrary electronic states. In this study, for brevity, in order to focus on a minimal description of the electron dynamics, we used the ground state electron density at an initial simulation time to construct a reference Hamiltonian for the electrons as the first stage of introducing our proposed scheme. During a simulation, Fock matrix was unchanged.
ii.2.2 Group localized orbitals of subgroup
First, we classify into subgroups, for example, monomers. Here, we can safely expect that the atom center of the ith Löwdin atomic orbital, , obtained by Löwdin orthogonalization is the same as that of the original ith AO orbital. In addition, arbitrary divisions of the component atoms among the total system are possible.
The block structure of the Fock matrix within the Löwdin basis set is expressed by
(5) 
where we used the property of as a realsymmetric matrix. Here, denotes the ith sub group while is the number of sub groups in the system.
The diagonalization of diagonal block corresponding to subgroup ,
(6) 
gives rise to the unitary transformation matrix, , of which the column vectors are the linear coefficient vectors of the localized eigenstates expanded in terms of the Löwdin orthogonalized atomic basis functions for the group G. The dagger symbol attached to a matrix indicates its adjoint form. is the diagonal matrix having eigen energies, , of the corresponding subgroup G as their elements, and is the number of local basis sets spanned at the group G. Here, . The collection of these localized group eigenstates over all groups, , spans the same Hilbert space of that of the original atomic orbitals.
In the equations above we implicitly indicated the complex Hermitian forms of the Fock operators and related matrices because, generally, they are complex and Hermitian. This is the case concerning spinorbit couplings. However, because we do not treat spinorbit couplings here, the Fock matrices in the AO and Löwdin orthogonalized AO basis representations are real and symmetric. Correspondingly, the transformation matrices for the construction of the group localized orbitals are real, orthogonal matrices.
The important feature is that these sets of localized orbitals are not orthogonalized between different groups; this provides a diabatic character in the representation with use of the collection of these orbital sets. Thus, we can obtain group localized orbital sets required for the GDF representation, as explained in the next subsection.
ii.2.3 Group diabatic Fock matrix
The group diabatic representation of the Fock operator is constructed as follows. The divided blocks in the Löwdin representation are transformed to the group interaction representation using ,
(7) 
which yields the GDF matrix,
(8) 
In this form, the submatrices in the diagonal blocks, , are the diagonal matrices corresponding to the local group eigen energies, while the submatrices placed at offdiagonal blocks, , describe the interactions with different blocks.
Note also that the Fock matrix in group diabatic representation does not lose any part of that in AO presentation, namely, no approximation is introduced during the transformation between these two representations.
ii.2.4 Transformation of observable, Fock and density matrices
A matrix representation of any observable operator in group diabatic basis set, , is related to that of the original AO basis set, , as follows:
(9) 
where
(10) 
The Fock matrix obeys the same transformation rule and is obtained by setting in the above equations, where we know and .
In contrast, the transformation of the density matrix between the original AO basis set and the group diabatic one is given by
(11) 
Note that the unitarity of assures total electron conservation for this transformation, .
ii.2.5 State couplings
In our theoretical method, the essential elements required for the construction of the lightelectron couplings are , , related to Eq. (9). Here, boldface denotes a vector in threedimensional Cartesian space, and denotes a composite variable of the electron position in threedimensional space. The first and second operators are responsible for the lightelectron couplings in length and velocity forms. FaisalLGVG (); ACPlaser () In this article, we neglect the nonadiabatic coupling and molecular motion to allow us to focus on the present electron dynamics scheme using the group diabatic representation.
ii.2.6 Time propagation of density matrix in GD representation
In the group diabatic representation explained in the previous subsection, the electron dynamics are naturally described in terms of the Liouvillevon Neumann equation and the associated density matrix as follows: where . Here, is a lightelectron coupling matrix. The lightelectron couplings are described by for the length gauge and for the velocity gauge. FaisalLGVG (); ACPlaser () Here, and are the threedimensional electric field vector and electromagnetic field vector potential, which generally depend on a point in a threedimensional space, respectively. Because the wavelength of light treated here is sufficiently large for the size of the molecular system treated in this article, we can safely employ the long wave length approximation that and do not depend on the spatial location and depend only on time. FaisalLGVG (); ACPlaser () In this study, we employed the length gauge. The time steps were set to 4 attosecond and we followed the dynamics for 15 femtosecond in all the calculations. We used the fourthorder Runge Kutta time integrator. If we directly consider the core orbitals, the spectral range of the GDF matrix becomes too large and a very small time step is required for a numerical convergence of total electron number as a trace of density matrix during a dynamics. As one way to avoid this, we first set to zero sufficiently small interactions between the core and other orbitals in the GDF matrix and by setting the original core orbital energies to a uniform value, such as the lowest valence orbital energy. This treatment is substantially equivalent to a neglection of core orbital dynamics and interaction with the valence orbitals. However, this does not reduce computational cost since the matrix dimension does not change. Of cause, we can extract the small size of valence and virtual space. In such a treatment, a computational cost will be reduced. But, this reduction treatment is not applied in the numerical demonstrations in the present article.
The electron matrices required for the dynamics calculation, namely, the Fock, dipole, and electron velocity moment matrices within the AO representation were obtained by using the electronic structure package, NTChem2013. NTChem ()
ii.2.7 Initial density matrix: local excitation and electron filling
The initial density matrices in the GDF representation are set so that diagonal elements in each diagonal block of the corresponding monomer should be occupied up to the number of electrons assigned to this monomer. The simple example will be shown later in this subsection. In a restricted case with the same spatial orbitals for different alpha and beta spins, the GDF local orbitals of each monomer are occupied by up to the half the number of electrons assigned to this group from the lowest energy orbital. Therefore, this initial density matrix, by construction, differs from that of the true ground state of the whole system and corresponds to the density matrix of a slightly excited state.
Based on this referential configuration in the group diabatic representation mentioned above, we can further introduce an excitation configuration measured from the referential occupations of the group diabatic localized orbitals. To prepare the initial densities in the cases of excess or deficient electrons in each monomer, we set the occupations to the corresponding number of electrons in each monomer.
For example, let us consider the spinrestricted case with three monomers each of which has 2 electrons and 3 orbitals. We suppose that we want to construct the initial condition such that the first and third monomers are initially excited from local HOMO to LUMO. In this case, the setting way of initial density matrix in GDF scheme is as follows: [1] set reference density matrix, namely GDF ground state, and [2] carry out local HOMO LUMO excitation and obtain the aimed GDF density, namely,
(12) 
For simplicity of presentations in this article, we employed the spinrestricted representation and the electron occupations of the neutral monomers at an initial simulation time. We also note that here we treat cases of overall singlet spin state in this article. The details and calculation results for the spinunrestricted and charged cases will be reported in our future articles.
By using this scheme, we can analyze the dynamical electron redistribution in molecular aggregate systems driven by the differences in the electron affinities between local groups and the opticalelectron interactions.
Iii Numerical applications
We examined two types of systems as a demonstration of the present scheme. The first example concerns the charge migration dynamics triggered by an external light field in a naphthalene(NPTL)tetracyanoethylene(TCNE) dimer, which used as a typical electron donoracceptor system; the results are used as the first illustration of our scheme. The second example is the unpaired electron dynamics in an excited 20mer ethylene system. This is used as a model of exciton transfer within a molecular ring system. Though the component monomers are different, such a ring shape is found in one unit of lightharvesting antenna systems LHanthenaring (); MOQME () where, an internal electron migration can affect an energy transfer among ring aggregates. The ring shape also found in a test system consisting of weakly interacting excited atoms at ultra cold temperatures as an example of exotic system for the experimental examination of the mechanism of quantum dynamics and its control by an external field associated with highly structured excited states. atomringRydspin (); Eisfeldultracoldcoins () In both cases, the time propagations of the density matrix were carried out in the Hilbert spaces spanned by all the AOs. In the present article, for simplicity, nonadiabatic couplings are omitted and cases of spin singlet state for the total systems are treated.
iii.1 Charge migration dynamics: NPTLTCNE dimer
Here, by applying the present theoretical method to a NPTLTCNE dimer system, we observed the initiation of charge separation triggered by the electronlight interaction, where the light had the ordinary wave length (700 nm) of sunlight and distinct polarization vectors along molecular axes. The balance of the photon energy and the energy gaps of localized orbitals of monomers are discussed in the supplementary material with respect to the mode analysis of electron dynamics. Because of the properties with respect to electron addition, as well as ejection (the adiabatic electron affinity of TCNE is large (3.16 eV) NIST (); adiaelecaffinityTCNE (), as is the adiabatic ionization potential of NPTL (8.14 eV) NIST (); ionizationenergyNPTL ()), the NPTLTCNE set makes an ideal donoracceptor molecular pair for checking the description of charge transfer with using the present GDF scheme.
In the Fig. 1, we present the geometry configuration of this dimer to show the relative orientations of monomers. Each monomer was optimized at the DFT/631G(d) level of theory using the CAMB3LYP exchange correlation functional.CAMB3LYP () The principal axis of NPTL is parallel to the Xaxis, while that of TCNE is parallel to the Yaxis. The molecular planes of these flat molecules are parallel to the XY plane. Although we treat a dimer here, we partly considered the reported crystal data crystalNTFLTCNE1967 () with respect to the relative orientation as explained in the figure caption. Note that the relative orientation employed here gives rise to a nonvanishing overlap between the frontier orbitals, that is, the HOMO of NPTL and the LUMO of TCNE, which is included in the panel (b) in Fig. 1.
Throughout this article we use a primitive combination of a basis set 631G(d) and functional CAMB3LYP because we focus on the validity of methodology. We note that that the dependency of the dynamics on the basis function and functional are important issue for charge transfer dynamics in a framework of DFT and TDDFT, which will be reported in our future articles.
By using this system, we examine (i) the charge migration from the initial neutral electron filling for monomers in the GD representation with and without a light field (ii) the dependency of the monomer distance on the dynamics with the different separations, (3.0, 3.5 and 4.0 Å) and (iii) the dependency of the field polarization directions on the results. We summarize the results of the charge dynamics of the two monomers in Fig. 2 and 3. The charges associated with the monomers were evaluated using Mulliken population analysisSzabo () by using a density matrix of the whole system. The charges of monomers were simply evaluated as a sum of those of constituent atoms in them. The distance between the molecular planes increases from the left column to the right column, while the laser field conditions vary with respect to the rows shown in the figure.
First we start from Fig. 2. The top three panels, (a/b/c), present the results in the cases without a light field. Note that, initially, each monomer is neutral within the GDF representation. From the timedependent charges of the monomers we found that there was a small amount of charge transfer oscillation between these two monomers. Note that negative and positive charges in the figure correspond to the acceptance and release of electrons, respectively. The amplitude is large when the distance is small, which is a reasonable result because the magnitude of molecular interaction depends qualitatively on the overlap of relevant orbitals.
Next, we see the cases where the field has a polarization vector along the Xaxis. These results are shown in the panels (d/e/f) in Fig. 2. In the panels (d) and (e), we find the sign of charge separation for the cases of distances of 3.0 and 3.5 Å , respectively, although not completely. As shown in panel (d), for the shortest separation of 3.0 Å, we found that the charge transfer was enhanced compared to the fieldfree case. The revival period after charge migration is about 8 fs. Interestingly, by increasing the separation to 3.5 Å, the charge separation signature appears as shown in the panel (e), where the half period is 14 fs. Note that the revivals are due to the fully coherent approach associated with a fixed molecular geometry. The observation here indicates that a reduction in the orbital overlap increases the time required to reach the peak of charge separation under an external field. In the other word, a very small separation is disadvantageous for meaningful charge separation between the donor and acceptor molecules, possibly arising from possible backandforward electron donation. However, as intuitively understood, a very large distance also reduces the efficiency of charge transfer because of the small interactions caused by the small overlaps of the relevant orbitals. As seen in the panel (f), which illustrates a case of a separation of 4.0 Å, the large distance suppresses the orbital overlap and, thus, the charge transfer time becomes too large. In fact, the charge transfer is negligible at this time scale.
In contrast, the optically forbidden direction of light polarization for the system does not induce the charge transfer. In fact, as shown in the panels (g/h/i) in Fig. 2, which correspond to the polarization direction along the Yaxis, we found negligible charge separation, similar to the fieldfree cases shown in (a/b/c) in Fig. 2.
Next, we see the results of the cases with a Z polarization of the light field which are shown in the panels (j/k/l) in Fig. 2. Comparing these results with the previous cases, (a–i) in Fig. 2, we found two clear different features in the charge migration dynamics. First, there is an increase in the number of oscillation in the envelope dynamics in (j) and (k), and, second, there is an enhancement in the moderate charge transfer at the largest distance, 4.0 Å as shown in the panel (l). As the distance increases from (j) to (l) the timedependent charge behavior becomes smooth.
Let us proceed to the cases with more than two components in the field polarization vectors, as shown in the panels (a–l) in Fig. 3. We will see that the charge separations are significant in the cases including both X and Zcomponents. First, the results in cases including both X and Ycomponents in the light polarization are presented in panels (a/b/c) in Fig. 3. Comparing with the results shown in (d/e/f) in Fig. 2, we find that the addition of the Ycomponent reduces the charge separation. In contrast, the combination of X and Z components drastically increase the efficiency of charge separation, as shown in panels (d/e/f) in Fig. 3. Among the cases presented here, this combination gave rise to the best performance. It may be instructive to compare with the case including only X component. According to the discussion in the cases of (d/e/f) in Fig. 2 with respect to the reduction of electron overlap with increase of the monomer distance, the difference in the degrees of charge separations between (k/l) in Fig. 2 including only z component in the poralization vectors and (e/f) in Fig. 3 including x and z components indicate that the increase in the revival period and the emergence of the associated charge separation can be attributed to the suppression of electron backdonation caused by the addition of the Xcomponent to the polarization vectors. Interestingly, the Ypolarization of light field in this molecular configuration suppressed the charge separation which can be read from the cases with the combination of Y and Zcomponents in the panels (g/h/i) in Fig. 3. Then, how about a charge migration dynamics between monomers under the use of the polarization vector including all the components? We displayed the results for the case including three components, X, Y and Z in the polarization vectors of light in the three panels (j), (k) and (l) in Fig. 3 with different distances of monomer planes, 3.0, 3.5 and 4.0 Å, respectively. These panels show the signature of charge separation induced by the external light field. Differently from the first two panels (j) and (k) with smaller distances of planes, in the panel (l) with the largest distance of 4.0 Å, we can see a clear charge separation in this time scale.
The understanding of the difference in charge separability with respect to the polarization vectors in the cases of largest distance, 4.0 Å corresponding to the panels (c/f/i/l) of Fig. 3 are supported by the information of the dynamics of induced dipole moment vectors and their Fourier transformation associated with the optical responsibilities of the system. The propensity of optical transition and its magnitude for the monomers are key factors. This discussion is given in the supplementary material including detailed explanations accompanied with the data of the transition dipole moment vectors and localized orbital energies of monomers. Generally speaking, in the present demonstrations, the reduction and increase in the charge separation observed here by changing these multi components in polarization vectors are attributed to a kind of electronicbath effect associated with superpositions of electronic mode oscillations having various time periods induced by simultaneous excitations of different states.
Through the investigation of the simple donoracceptor system, we obtained the following two findings; (I) a timedependent approach coupled with the group diabatic representation yields a microscopic information on electron properties such as charge dynamics, which can not be obtained by a static analysis or nuclear wave packet calculation ( The statement above does never deny static and nuclear wave packet approaches. In fact, they can offer a useful and high level information on electronic properties of charge transfer and nuclear quantum effect. ) , and (II) the initiation dynamics of the charge separation are sensitively dependent on the distances of the constituent monomers and the optical properties associated with the donoracceptor system, which also cannot be extracted from a static analysis.
Of course, the effects of electronic nonadiabaticity, as well as molecular dynamics, are also important for the charge separation dynamics. However, here, we have not touched these issues for clarity. They will be discussed and reported in our future articles.
iii.2 Unpaired electron dynamics: 20mer ethylene
In the next demonstration, we focus on the unpaired electrons in excited electron dynamics. The diffusion dynamics of the excitons over weakly interacting monomers are examined. We examine the differences between results obtained with and without the group localization procedure. The group localization scheme provides us with a clear view of the migration dynamics of local excitons in the monomers in the presence and absence of a laser field.
iii.2.1 Effective unpaired electron (EUPE)
First, we define an unpaired electron as used in this paper. An effective unpaired electron density matrix was constructed by extracting the components corresponding to an occupation of almost one (half filling) of the natural orbitals from the original density: unpaired (); nagashima2012 ()
(13) 
where with and being the ith natural orbital and its corresponding natural population, respectively. Here, denotes a oneelectron density operator for the whole system. This treatment enables us to obtain an information concerning the polyradical features of the complex system. b12nad ()
For example, in the GD representation case, quantities of unpaired electrons for monomers are evaluated as follows:

Diagonalize electron density matrix, , in the group diabatic representation, where and are a unitary matrix and diagonal natural population matrix.

Create a unpaired electron density matrix,

Convert into the form in AO representation and perform the Mulliken population analysis for monomers.
iii.2.2 System information on the 20mer ethylene
For the assessment of the present method, we treated an aggregated system consisting of a circle of twenty ethylene molecules, as shown in Fig. 4. The centers of mass of the monomers are uniformly placed on the circle, which has a radius being 12 Å. The CC lines of planer ethylene molecules are aligned vertically to the circle plane. The normal vector of each ethylene is parallel to the tangent vector of this circle at the position of the monomer.
The geometrical structure of ethylene was determined at the CAMB3LYP/631G(d) level of theory. We also employed the same level of abinitio calculation for the 20mer ethylene system as that used for ethylene monomer.
iii.2.3 Exciton migration dynamics
Here, we show that the group diabatic representation is an advantageous for investigating the exciton migration dynamics in molecular aggregate systems containing weak interactions between monomers.
First, we present the results in the cases without any group diabatization. Note that, here, we employ in the GDF representation. This means that the 20mer ethylene system was treated as one monomer system. Initially, the entire system was excited from the HOMO to the LUMO and the corresponding electron density matrix was prepared. We note that the treatment of the total system as a monomer is rather natural. This first presentation treating a whole system as a monomer is aimed to compare with the local excitation analysis by using group diabatic scheme performed later in this article. This means to see how the introduced method is useful in the local diabatic analysis in a time domain, which is the main topic in this article. The panels (ac) of Fig. 5 summarizes the results for the cases without any group diabatization.
For this initial electronic density matrix, two types of dynamics simulations were carried out: One with and one without an external optical field.
In the cases with light fields, as model cases of modification of excited electron transfer stirred by solar light, we employed the same parameters for the continuum light field except for the polarization vectors, which were used as those in the case of the NPTLTCNE dimer. The details of parameters are included in the figure captions. The polarization direction vectors of the light field used here have inplane and outofplane components with respect to the plane on which the center of masses of the ethylene monomers are positioned.
Note also that the symmetry of the system is broken in a case accompanied with this light field, for example, symmetry operations with respect to the XZ plane and this will also be true for the induced exciton migration.
The number of unpaired electrons in each monomer was evaluated via Mulliken population analysis for the unpaired electron density, Eq. (13). For clarity, we have used the label number of monomer as the origin in plotting the timedependent data sequences.
As shown in the panel (a), in the absence of a light field, the unpaired electrons are distributed uniformly over the twenty monomers and no migration occurs. This indicates that the unpaired electron affinities are balanced among these moieties in this setting of calculation. In contrast, unpaired electron dynamics between the monomers appear weakly in the presence of an external radiation field, as shown in the panel (b) and (c). The field polarization vectors in these two cases, are related by inversion with respect to the XZ plane. Thus, the patterns of unpaired electron dynamics is in an inverse relationship between these two cases for the XZ plane.
From the resultant uniform excitation of the monomers shown in (a–c), the initial excitation of the entire system is not so useful for the analysis of the time dependent monomer interactions. The difficulty in the interaction analysis without group diabatization is resolved by applying a group diabatized representation and the associated local initial excitation explained in the Sect. II. Now, let us proceed to trial calculations using the group diabatic representation discussed in this article.
Next, we examine the case where initially only one ethylene molecule is excited from the HOMO to the LUMO in the localized canonical orbitals of this moiety. Here, group localizations are carried out by using all the monomers. This means that we employed in the GDF representations. The results are presented in the panels (d)–(f) of Fig. 5.
As seen in the panel (d), in the absence of an external radiation field, the highly localized unpaired electrons in the 1st monomer are separately transferred to 11th monomer at 10 fs, accompanied by a moderate broadening in the distribution over the monomers. Symmetry with respect to the XZ plane was also observed.
In the case with the external light field of the panel (e), we found that the exciton transfer becomes slow associated with the field induced change in the interactions between monomers. The symmetry of the migration dynamics of the unpaired electrons with respect to the XZ plane is clearly broken by the light polarization direction. We also present the results, in the panel (f), of the calculation with a symmetric light polarization compared to the case of (e) with respect to the XZreflection plane. In this case of (f), the symmetry breaking is the inverse of that (e).
As a final example of unpaired electron dynamics, we examine the cases of initial twosite single excitations from the local HOMO to the LUMO for the 1st and 11th monomers, for which results are included in the panels, (g), (h) and (i). The conditions of the light fields used in the panels (g/h/i) are the same as those in (d/e/f) in this ordering, for which details are included in the figure caption. As shown in panel (g), there is a clear bifurcation and confluence of two localized excitions accompanied with moderate dispersion. These dynamics can be modified by an external light field. The panels, (h) and (i), show the two cases with polarized directions related by reflection with respect to the XZ plane. From these data, we can observe that the directions of rectification of the exciton flows induced by the light field are inverse with respect to the XZ plane for these two cases. The difference in the Fourier spectrum and timedependent induced dipole moment vectors between one and two initial local excitation cases are summarized in the supplementary material, where the readers can find the information on the mode of electronic motions dominantly triggered by the initial local excitation and following optical response of electrons. As seen in the demonstration of initial twoFrenkelexciton case, the method can safely describe the details of the exciton migration dynamics over the monomers. This is also the case for the one exciton case as shown above. We give schematic movies in the supplemental materials for the two excitation cases with and without the light field, where readers can see the unpaired electron wave packet dynamics, especially the extent of the delocalization of packets, time scale in it and rectification feature associated with external light field.
So far we examined unpaired electron dynamics. How about the dynamics of charge of local sites ? In Fig. 6, we showed the result of Mulliken charge migration dynamics over twenty monomers. The panels of (a)(c) display the result of the dynamics starting from initial density matrix corresponding to an excitation of 1st monomer from local HOMO to LUMO while in the panels of (d)(f) 1st and 11th monomers are initially excited in the same way. Here the diabatic group number, , is twenty. The field parameters of (a)(f) in Fig. 6 are the same as those of (d)(i) in Fig. 5 in this order. The plotting manner are the same as that of Fig. 5.
By comparing the charge migration dynamics in the panels (a)(f) in Fig. 6, with the time dependent unpaired electrons of local sites shown in the corresponding panels (d)(i) in Fig. 5, we find the following features and obtain the results from them:

As a summary, the charge migration clearly seen in Fig. 6 indicates that the initially prepared local excitons are characterized by a charge transfer type of exciton. The charge density wave packets slightly delocalized over two or three monomers propagate spatially along the molecular chains in the molecular aggregate system, and drive the charge migration over twenty ethylene monomers. We can consider the present electron migration dynamics as a charge transfer excition type of dynamics.

In fact, as seen in the panel (a) of Fig. 6 with the one site excitation and without the light field, the initial Mulliken charge of the initially excited 1st ethylene site is zero and this monomer, soon after the onset of the dynamics, obtains the positive charge from the nearest two 2nd and 20th sites. Then this positive charge stored in 1st site in turn propagates in a splitting manner via 2nd and 20th monomers toward the 11th monomer at the counter position with respect to the 1st monomer. This two charge density wave packets encounter at 11th around the time of 12 fs accompanied with a constructive interference.

The panel (d) of Fig. 6 for the case of initially prepared two excitations show the same feature as (a) of the same figure except for the double numbers of spawned charge density wave packets.

The packets have the width over two or three sites during dynamics, which characterize the exciton delocalization. This means that the initially prepared Frenkel like exciton at each local site was converted into stable delocalized excitons accompanied with a charge during the dynamics.
Though these calculations include rather artificial initial setting of local excitation, we can extract inherent dynamical trend in this molecular aggregated system between local excitation and local charge by using GDF electron dynamics method.
Thus, through the dynamics calculations with and without the group diabatization representation, various local excitations and laser fields, we found that this analysis has a utility to obtain an information on a dynamical propensity of excited electrons in aggregated systems having sparse networks of the interactions between the constituent monomers.
Iv Summary
In this paper, we have introduced a calculation scheme for excited electron dynamics based on the group diabatic Fock representation. We verified that this GDF electron dynamics method allows for the concise description and analysis of the excited electron migration dynamics in molecular aggregate systems. This was assessed by using an elemental light energy conversion material made of electron donor and acceptor molecules and a onedimensional system consisting of circularly oriented monomers. The dynamics were characterized according to the inherent gradient of electron affinities among the local molecular groups under the employed conditions of initial excitations and external laserelectron couplings.
The present scheme is advantageous for the future ab initio modeling of excited electron migration dynamics in a moderately large system. This is because the GDF representation provides a clear strategy for the extraction of the active orbitals at each local group site by setting the orbital energy range related to the excited electron transfer under consideration.
In this aspect, the promising scheme recently developed by Shimazaki, Kitaura, Fedrov and Nakajima as a fragment type duallayer self consistent field theory for the treatment of a large sparse system GMO () will play an important role in the exploration of the roles of the structured but complex interactions of many types of molecular aggregate systems. In fact, the present work is partly inspired by the work mentioned above and can be combined with it by replacing the localization scheme with e.g. BoysBoys (), PipekMezeyPM () methods and so on.
The external, as well as internal, fields imposed on the systems affect the excited electron dynamics and energy transfer. Therefore, the molecular motion associated with nonadiabatic transition among complex excited states are also important at longer timescales than that considered in the present article. envassistqweet () Furthermore, in systems containing metal atoms, the spinorbit couplings become important for excited state dynamics involved with intersystem crossings. DanielACRSOCQMDYN (); FedrovAIMSSOC () We will report these issues in future articles by using the extended version of the present GDF electron dynamics scheme.
V Supplementary Material
See supplementary material for Fourier analysis of induced dipole moment of the NPTLTCNE dimer and 20mer ethylene systems, the fundamental properties of transition dipole moments as well as localized orbital energies for monomers and the schematic movies of dynamics of local unpaired electrons in the latter system.
Acknowledgements.
The authors thank Dr. Michio Katouda and Dr. Keisuke Sawada for their advises, discussions and technical support concerning the use of the NTChem code and parallel calculations. We are grateful to Prof. Kazuo Kitaura and Prof. Unpei Nagashima for the valuable information on the molecular interaction concept in terms of molecular orbitals. We also appreciate the discussions with Dr. Tomomi Shimazaki, Prof. Yuzuru Imamura, Prof. Motomichi Tashiro and Prof. Mikiya Fujii regarding the charge separation mechanism in solar cell systems. This research was supported by MEXT, Japan, “NextGeneration Supercomputer project” (the K computer project) and “Priority Issue on PostK computer” (Development of new fundamental technologies for highefficiency energy creation, conversion/storage and use). The part of computations in the present study was performed using Research Center for Computational Science, Okazaki, Japan.References
 (1) V. May and O. Khn, Charge and Energy Transfer Dynamics in Molecular Systems, (WileyVCH, Berlin, 2000).
 (2) G. D. Scholes and G. Rumbles. Nat. Mater. 5, 683 (2006).
 (3) S. H.Schiffer and A. A. Stuchebrukhov, Chem. Rev. 110, 6939 (2010).
 (4) N. J. Hestand, R. V. Kazantsev, A. S. Weingarten, L. C. Palmer, S. I. Stupp, and F. C. Spano, J. Am. Chem. Soc. 138, 11762 (2016).
 (5) M. R. Wasielewski, Acc. Chem. Res. 42, 1910 (2009).
 (6) C. J. Brabec, G. Zerza, G. Cerullo, S. De Silvestri, S. Luzzati, J. C. Hummelen, and S. Sariciftci, Chem. Phys. Lett. 340, 232 (2001).
 (7) R. Huber, J.E. Moser, M. Grätzel, and J. Wachtveitl, J. Phys. Chem. B 106, 6494 (2002).
 (8) M. Mohseni, P. Rebentrost, S. Lloyd, and A. AspuruGuzik, J. Chem. Phys. 129, 174106 (2008); P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. AspuruGuzik, New J. Phys. 11, 033003 (2009).
 (9) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
 (10) M. A. L. Margues, A. Castro, G. F. Bertsch, and A. Rubio, Comput. Phys. Comm. 151, 60 (2003).
 (11) M. Noda, K. Ishimura, and K. Nobusada J. Comput. Phys. 265, 145 (2014).
 (12) K. Takatsuka and T. Yonehara, Phys. Chem. Chem. Phys. 13, 4987 (2011).
 (13) T. Yonehara, K. Hanasaki, and K. Takatsuka, Chem. Rev. 112, 499 (2012).
 (14) K. Takatsuka, T. Yonehara, K. Hanasaki, and Y. Arasaki, Chemical Theory beyond the BornOppenheimer Paradigm (World Scientific, 2015).
 (15) Y. Xie, J. Zheng, and Z. Lan, J. Chem. Phys. 142, 084706 (2015).
 (16) A. Cembran, L. Song, Y. Mo, and J. Gao, J. Chem. Theory Comput. 5, 2702 (2009).
 (17) B. Shi, F. Gao, and W. Liang, Chem. Phys. 394, 56 (2012).
 (18) J. Li, I. Kondov, H. Wang, and M. Thoss, J. Phys. Chem. C 114, 18481 (2010).
 (19) A. Szabo and N. S. Ostlund, Modern quantum chemistry: Introduction to Advanced Electronic Structure Theory (Dover, New York, 1996).
 (20) T. Yonehara and K. Takatsuka, J. Chem. Phys. 132, 244102 (2010).
 (21) T. Yonehara and K. Takatsuka, J. Phys. Chem. A 117, 8599 (2013).
 (22) Y. Suzuki and K. Watanabe, Phys. Rev. A 94, 032517 (2016).
 (23) A. Abedi, N. T. Maitra, and E. K. U. Gross, Phys. Rev. Lett. 105, 123002 (2010).
 (24) P. Krause, T. Klamroth, and P. Saalfrank, J. Chem. Phys. 127, 034107 (2007).
 (25) M. Amano and K. Takatsuka, J. Chem. Phys. 122, 084113 (2005).
 (26) J. Zanghellini, M. Kitzler, C. Fabian, T. Brabec, and A. Scrinzi, Laser Phys. 13, 1064 (2003).
 (27) T. Kato and H. Kono, Chem. Phys. Lett. 392, 533 (2004).
 (28) T. Kunisada, H. Ushiyama, and K. Yamashita, Chem. Phys. Lett. 653, 345 (2015).
 (29) H. D. Meyer, F. Gatti, and G. A. Worth, eds., Multidimensional Quantum Dynamics: MCTDH Theory and Applications, (WileyVCH, Weinheim, 2009).
 (30) T. Sato and K. K. Ishikawa, Phys. Rev. A 88, 023402 (2013).
 (31) H. Miyagi. L. B. Madsen, Phys. Rev. A 87, 062511 (2013); ibid. 89, 063416 (2014).
 (32) I. S. Ulusoy and M. Nest, J. Am. Soc. Chem. 133, 20230 (2011).
 (33) U. Schollwöck, J. Phys. Soc. Jpn. 74, 246 (2005); A. J. Daley, C.Kollath, U. Schollwöck and G. Vidal, J. Stat. Mech. 4, 04005 (2004); E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. 101, 140601 (2008).
 (34) T. Yonehara and K. Takatsuka, Chem. Phys. 366, 115 (2009).
 (35) H. Ushiyama and K. Takatsuka, Angew. Chem. Int. Ed. 46, 587 (2007).
 (36) K. Nagashima and K. Takatsuka, J. Phys. Chem. A 116, 11167 (2012).
 (37) K. Yamamoto and K. Taktatsuka Chem. Phys. 475, 39 (2016).
 (38) T. Yonehara and K. Takatsuka, J. Chem. Phys. 144, 164304 (2016).
 (39) X. Li, D. J. Haxton, M. B. Gaarde, K. J. Schafer, and C. W. McCurdy, Phys. Rev. A 93, 023401 (2016).
 (40) T. Matsuoka and K. Taktatsuka, J. Chem. Phys. 146, 134114 (2017).
 (41) K. C. Kulander, Phys. Rev. A 36, 2726 (1987).
 (42) S. Meng and E. Kaxiras, J. Chem. Phys. 129, 054110 (2008).
 (43) A. Castro, M. A. L. Marques, and A. Rubio, J. Chem. Phys. 121, 3425 (2004).
 (44) C.Lun, J. S. Evans, and T. V. Voorhis, Phys. Rev. B 74, 155112 (2006).
 (45) Fundamentals of TimeDependent Density Functional Theory M. A. Marques, N. T. Maitra, F. M. Nogueira, E. K. U. Gross, and A. Rubio, (Eds.) ( Springer, 2012 ).
 (46) T. Minami, H. Fukui, H. Nagai, K. Yoneda, R. Kishi, H. Takahashi, and M. Nakano, J. Phys. Chem. C 113, 3332 (2009); R. Kishi, M. Nakano, T. Minami, H. Fukui, H. Nagai, K. Yoneda, and H. Takahashi, J. Phys. Chem. A 113, 5455 (2009); M. Nakano, R. Kishi, T. Minami, and K. Yoneda, Molecules 14, 3700 (2009).
 (47) H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003).
 (48) F. H. M. Faisal, Phys. Rev. A 75, 063412 (2007).
 (49) K. Takatsuka and T. Yonehara, Adv. Chem. Phys. 144, 93 (2010).
 (50) T. Nakajima, M. Katouda, M. Kamiya, and Y. Nakatsuka, Int. J. Quantum Chem. 115, 349 (2015).
 (51) T. Mirkovic, E. E. Ostroumov, J. M. Anna, R. van Grondelle, Govindjee, and G. D. Scholes, Chem. Rev. 117, 249 (2017).
 (52) H. Labuhn, D. Barredo, S. Ravets, S. de Léséleuc, T. Macri, T. Lahaye, and A. Browaeys, Nature 534, 667 (2016).
 (53) S. Wüster, A. Eisfeld, and J. M. Rost, Phys. Rev. Lett. 106 153002 (2011).
 (54) http://webbook.nist.gov/
 (55) D. Khuseynov, M. T. Fontana, and A. Sanov, Chem. Phys. Lett. 550, 15 (2012).
 (56) M. C. R. Cockett, H. Ozeki, K. Okuyama, and K. Kimura, J. Chem. Phys. 98, 7763 (1993).
 (57) T. Yanai, D. P. Tew, and N. C. Handy, Chem. Phys. Lett. 393, 51 (2004).
 (58) R. M. Williams and S. C. Wallwork, Acta. Cryst. 22, 899 (1967).
 (59) K. Takatsuka, K. Yamaguchi, and T. Fueno, Theor. Chim. Acta 48, 175 (1978).
 (60) T. Shimazaki, K. Kitaura, D. G. Fedorov, and T. Nakajima, J. Chem. Phys. 146, 084109 (2017).
 (61) S. F. Boys, Rev. Mod. Phys. 32, 296 (1960).
 (62) J. Pipek and P. G. Mezey, J. Chem. Phys. 90, 4916 (1989).
 (63) J. Eng, C. Gourlaouen, E. Gindensperger, and C. Daniel, Acc. Chem. Res. 48, 809 (2015).
 (64) D. A. Fedorov, S. R. Pruitt, K. Keipert, M. S. Gordon, and S. A. Varganov, J Phys. Chem. A 120, 2911 (2016).
Supplementary material for “A quantum dynamics method for excited electrons in molecular aggregate system using a group diabatic Fock matrix”
Appendix A Suppl.1: Mode analysis of electron dynamics in NPTLTCNE dimer
In Fig. 7, we summarized the mode analysis of GDF electron dynamics in the dimer NPTLTCNE system in the main text of the article for the cases including the light field having 700nm wave length corresponding to one photon energy 0.065 Hartree, and normalized polarization vectors of (x,y,z) = , , and . Note that this light field is rather nonresonant to the typical local electronic transitions because the energy gaps of local HOMOLUMO being 0.0010(0.2721)=0.2731 for NPTL and 0.1278(0.3815)=0.2537 for TCNE are far larger than 0.065. The results in the cases with the distance of molecular planes being 4.0 Å are selectively shown, because they include key information of the response dynamics of electrons with respect to the light polarization vectors having the elements more than two in this molecular orientation against for the laboratory frame in the present demonstration.
Here, the time sequences of induced dipole moment and the Fourier transformation of them are displayed. For the light field, only the case of the polarization vector, , is shown for the guide of eye to show the frequency and time dependent behavior of the external radiation fields. An induced dipole of total system used for analysis was defined as the difference of the time dependent dipole moment measured from the dipole moment at the initial simulation time, with running over x, y and z. The broadenings of the peaks in FT spectrums are caused by the finite resolution associated with the simulation time, 15 fs.
The height for each component in FT spectrum at the peak of 0.065 a.u. in the frequency domain shown in the left column has an information on the extent of the response of electrons for the polarization component in light field. The panels tell us how the electrons respond to the light field, with respect to the induced dipole moment.
We can see that the light fields having the polarization vector including z component cause the gradual increase of z component in dipole moment as shown in the panels of (g) and (i) except for the panel (h) of the yz pair polarization. The low frequency peaks far below 0.065 in the panels (b) and (d) correspond to the slow increase of z component in the dipole moment associated with the unidirectional electron migration between monomers along z axis as shown in the panels of (g) and (i), respectively. We can see again that z component in the light field is essential for the electron transfer. In fact, there is no growth of dipole moment in the z direction in the panels (a) and (f) for the case including x and y components in the polarization vector but not having z component in it.
Interestingly, y component in light polarization vector suppresses the electron transfer along z axis caused by z component of the light polarization vector. As found in the panels (c) and (h) including y and z components in the polarization vector, the response of electrons along y axis is clearly superior to that along z axis. This indicates that for this field the system has a strong electronic transition property along y axis, which resultantly prevents the opportunities of the electron transfer along z axis parallel to the normal vectors of two molecular planes. This picture is supported by the fact that for the NPTL monomer as the electron donner the local HOMOLUMO transition dipole moment vector has large values for the y component and zero values for x and z component as seen in Tab. 1.
Appendix B Suppl.2:Mode analysis of electron dynamics in 20mer circle of ethylene molecule
Fig. 8 gives the electronic mode analysis in the cases of a 20mer circle of ethylene molecule with initial one and two excitations with and without light field. In the case of light, we employed the continuum light having 700nm wavelength (corresponding to one photon energy, 0.065 a.u.) and normalized polarization vector of (x,y,z) = . See also the main article. The induceddipole moment vectors and its Fourier transformation are displayed in the same manner as the previous section in this supplementary material. We find following features which characterize the present dynamics:

Excitation of electronic motion by the light of z polarization is significantly large compared to those of x and y polarization, which is consistent with the nonvanishing z component in the transition dipole moment vector with respect to the local HOMO and LUMO of the representative monomer, 1st ethylene, as shown in Tab. 3.

Electronic motion of z component experience multi photon excitation, which is observed from the spectrum peak around 0.4 a.u. in panels (c) and (d) and corresponding high oscillation behavior of z component of the dipole moment in the time domain as found in the panels of (h) and (i). This spectrum peak around 0.4 a.u. is consistent with the local HOMOLUMO energy difference, 0.0907(0.3140) = 0.4047 Hartree, of which values are shown in the caption of Tab. 3. And also, 700 nm corresponding to the one photon energy 0.065 Hartree is substantially nonresonant to the this typical electron transition.

The case of one local excitation is characterized by the deviation of x and y component as seen in the panels (h). This is reflected in the low frequency spectrum at the position almost near to zero (far lower than 0.065) as seen in the panel (c). The statement that this low frequency mode is attributed to the exciton migration is assisted by the fact that in the lightfield free case starting from the initial one site excitation we can find the low frequency mode in the spectrum around zero in the panel (a) and the slow dynamics of the induced dipole moment in the corresponding panel (f).

The low frequency peak seen in the panel (a) does not appear in the case of two local initial excitations as shown in the panel (b). This is cancellation effect due to the symmetry of the two unpaired electron wave packets starting from the 1st and 11th sites at the counter position in the circle. Note that here the dipole moments of total system is analyzed for simplicity. See also the schematic movie of this dynamics, of which explanation is given in the supplemental material. This cancellation effect remains to some extent under the shine of light field as seen in the panels of (d) and (i) though the symmetry is weakly broken as seen in the figure in the main article and the movie.
Appendix C Suppl.3:Schematics movies of unpaired electron dynamics in 20mer circle of ethylene molecule
Here we provide the schematic movies of dynamics of unpaired electron in the system of 20mer ethylene system, as mov20mertwoexcitons1.gif and mov20mertwoexcitons2.gif. These two files correspond respectively to the cases of initial two local site excitations without and with light field having the strength being 0.02, wave length of 700nm, and polarization vector of . The details for plotting data are explained below.
We constructed the movies by gathering the time sequence of the pictures in which the following two dimensional time dependent function, , corresponding to the quantity of unpaired electrons on sites are plotted in the xy plane on which all the center of masses (COMs) of ethylene monomers are placed,
(14) 
where
(15) 
the positions of COM of ith monomer,
(16) 
and being the quantity of unpaired electrons in ith monomer at a time t. Here, we set and (Å).
Appendix D Suppl.4:Transition dipole moment vectors and group localized orbital energies of monomers in the group diabatic representation
The transition dipole moment vectors are presented for NPTL and TCNE in the dimer system and 1st ethylene in the 20mer circle of ethylene molecules respectively in Tab. 1, 2 and 3. The associated local orbital energies are also given from the local HOMO2 to LUMO+2 in the captions of these tables. Note that the local properties of monomers presented here are the part of the group blocks in the matrix or vectors in the GD representation. With respect to the HOMO and LUMO energies of monomers, we also included in the captions the values obtained in the calculations for the totally isolated systems.
x  32  33  34  35  36  37 

32  1.0435  0.0000  2.5248  0.0000  0.0418  0.0000 
33  0.0000  1.2461  0.0000  2.1075  0.0000  0.0478 
34  2.5248  0.0000  1.2746  0.0000  2.0856  0.0000 
35  0.0000  2.1075  0.0000  1.0585  0.0000  2.5170 
36  0.0418  0.0000  2.0856  0.0000  1.1197  0.0000 
37  0.0000  0.0478  0.0000  2.5170  0.0000  1.3105 
y  32  33  34  35  36  37 
32  0.0000  0.0040  0.0000  0.0736  0.0000  1.3900 
33  0.0040  0.0000  0.0328  0.0000  1.4397  0.0000 
34  0.0000  0.0328  0.0000  1.4510  0.0000  0.0735 
35  0.0736  0.0000  1.4510  0.0000  0.0231  0.0000 
36  0.0000  1.4397  0.0000  0.0231  0.0000  0.0346 
37  1.3900  0.0000  0.0735  0.0000  0.0346  0.0000 
z  32  33  34  35  36  37 
32  3.7814  0.0000  0.0004  0.0000  0.0010  0.0000 
33  0.0000  3.7795  0.0000  0.0002  0.0000  0.0006 
34  0.0004  0.0000  3.7801  0.0000  0.0011  0.0000 
35  0.0000  0.0002  0.0000  3.7773  0.0000  0.0036 
36  0.0010  0.0000  0.0011  0.0000  3.7803  0.0000 
37  0.0000  0.0006  0.0000  0.0036  0.0000  3.7759 
x  30  31  32  33  34  35 

30  1.0762  0.0263  0.0000  1.6083  0.0000  0.0036 
31  0.0263  1.0362  0.0000  0.0052  0.0000  1.3513 
32  0.0000  0.0000  1.1384  0.0000  0.0025  0.0000 
33  1.6083  0.0052  0.0000  1.1767  0.0000  0.0035 
34  0.0000  0.0000  0.0025  0.0000  1.2254  0.0000 
35  0.0036  1.3513  0.0000  0.0035  0.0000  1.1979 
y  30  31  32  33  34  35 
30  0.0000  0.0000  0.0546  0.0000  0.0033  0.0000 
31  0.0000  0.0000  0.0048  0.0000  0.0107  0.0000 
32  0.0546  0.0048  0.0000  2.1112  0.0000  0.0039 
33  0.0000  0.0000  2.1112  0.0000  0.0104  0.0000 
34  0.0033  0.0107  0.0000  0.0104  0.0000  2.7252 
35  0.0000  0.0000  0.0039  0.0000  2.7252  0.0000 
z  30  31  32  33  34  35 
30  3.7800  0.0504  0.0000  0.0004  0.0000  0.0006 
31  0.0504  3.7813  0.0000  0.0005  0.0000  0.0000 
32  0.0000  0.0000  3.7827  0.0000  0.1305  0.0000 
33  0.0004  0.0005  0.0000  3.7868  0.0000  0.0451 
34  0.0000  0.0000  0.1305  0.0000  3.7833  0.0000 
35  0.0006  0.0000  0.0000  0.0451  0.0000  3.7815 
x  6  7  8  9  10  11 

6  22.6769  0.0000  0.0000  0.0000  0.0211  0.8328 
7  0.0000  22.6762  0.0000  0.0000  0.0000  0.0000 
8  0.0000  0.0000  22.6798  0.0000  0.0000  0.0000 
9  0.0000  0.0000  0.0000  22.6812  0.0000  0.0000 
10  0.0211  0.0000  0.0000  0.0000  22.7776  1.8637 
11  0.8328  0.0000  0.0000  0.0000  1.8637  22.5812 
y  6  7  8  9  10  11 
6  0.0000  0.0000  0.0430  0.0000  0.0000  0.0000 
7  0.0000  0.0000  0.0000  0.0031  0.0000  0.0000 
8  0.0430  0.0000  0.0000  0.0000  0.0895  0.0058 
9  0.0000  0.0031  0.0000  0.0000  0.0000  0.0000 
10  0.0000  0.0000  0.0895  0.0000  0.0000  0.0000 
11  0.0000  0.0000  0.0058  0.0000  0.0000  0.0000 
z  6  7  8  9  10  11 
6  0.0000  0.0003  0.0000  0.0000  0.0000  0.0000 
7  0.0003  0.0000  0.0000  0.0000  0.0254  1.0482 
8  0.0000  0.0000  0.0000  1.3581  0.0000  0.0000 
9  0.0000  0.0000  1.3581  0.0000  0.0000  0.0000 
10  0.0000  0.0254  0.0000  0.0000  0.0000  0.0000 
11  0.0000  1.0482  0.0000  0.0000  0.0000  0.0000 