A quantum dynamics method for excited electrons in molecular aggregate system using a group diabatic Fock matrix

A quantum dynamics method for excited electrons in molecular aggregate system using a group diabatic Fock matrix

Takehiro Yonehara takehiro.yonehara@riken.jp    Takehito Nakajima nakajima@riken.jp RIKEN Advanced Institute for Computational Science, Kobe 650-0047, Japan
July 25, 2019

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 time-dependent excitation of local sites in complex systems. In addition, light-electron couplings are considered. The present scheme is intended for investigations on the migration dynamics of excited electrons in light-induced energy transfer systems. The scheme was applied to two systems: a naphthalene(NPTL)-tetracyanoethylene(TCNE) dimer and a 20-mer 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. Kuhn-CET (); Scholes-rev-excitons (); Schiffer-Rev-PCET (); Spano-EET () 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 hole-electron 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. art-ps (); rt-exp-pcbm (); rt-exp-Graetzel (); env-assist-qw-eet ()

Theories for the treatment of real-time electron dynamics in complex excited states accompanied by the changes in chemical bonding, have been developed and applied, for example, real-time time dependent density functional theory, time dependent configuration interaction, and configuration state function schemes. RG-TDDFT (); octpus-rt-tddft (); gceed (); PCCP-eldt-KT-TY (); Nonadiabaticity-ChemRev (); Chemical-theory-beyond-BO-paradigm () These theories provide an excellent description of electron dynamics in total systems. However, there is no explicit scheme for the analysis of the time-dependent 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 time-consuming, 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 light-matter interactions.

This treatment enables us to treat arbitrary types of initial local excitation in a complex system. We have used the well-known group diabatic matrix representation of Fock operator. Lan (); Gao (); Shi (); Thoss-TiO2 () 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 light-matter 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 light-energy to chemical functionality via molecular aggregate systems. We show the efficiency of the introduced theoretical method by two example applications, a typical donor-acceptor 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, PSANB-laser (); en-coins-JPCA (); Bohmian-Suzuki-2016 (); Factorization-EN-Gross () 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 multi-electron wavefunction as a time-dependent linear combination of configuration state function constructed from a molecular orbital set and (b) a time-dependent 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 multi-state multi-electron wavefunctions

The time-dependent configuration interaction theory (TDCI) including the time dependent configuration state function (TDCSF) Saalfrank-JCP2007-TDCI (); TDCSF-Amano (); Nonadiabaticity-ChemRev (); Chemical-theory-beyond-BO-paradigm () and multi configurational time-dependent Hartree-Fock (MCTDHF) schemes mctdhf-scrinzi (); mctdhf-kato-kono () 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 multi-electronic 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 time-dependent 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. TDCSF-Amano (); Nonadiabaticity-ChemRev (); TDCSF-Kunisada-Ushiyama ()

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, book-Meyer-MCTDH () both the configuration coefficients and orbital parameters are simultaneously propagated according to the equation of motion derived from the application of the time-dependent variational principle to the quantum action. Here, time-dependent 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 one-electron self-consistent field, in MCTDHF-type of schemes, we can save the dimensions of the multi-electron configuration functions. This approach is usually applied to high-accuracy studies of field induced electron dynamics in small molecular systems with high accuracy. As practical schemes accompanied by the modelling of a time-dependent active orbital space in a sophisticated manner, the time-dependent complete active space self-consistent field (TDCASSCF) tdcascf-sato () theory and the restricted active space counterpart (TDRASSCF) tdrasscf-Madsen () 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 multi-electron adiabatic or pseudo-diabatic 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 post-Hartree-Fock theories. TD-CAS-CISD-Nest (); TD-DMRG ()

Applications with using (a) include a field-induced electron dynamics at the attosecond timescale, Diborane-laser (); PCCP-eldt-KT-TY () nonadiabatic electron wave packets in a proton-coupled electron transfer processes, Ushiyama-ACIE2006-PCET (); nagashima2012 (); TDCSF-pcet-yama-taka () the characterization of electron dynamics in highly quasi-degenerate excited states b12-nad (), and ultrafast photoionization process. MuCurdy-MCTDHF-appl-ionization (); TDCSF-ionization-matsu-taka () 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 Hartree-Fock (TDHF) and, Kulander-PRA-rt-TDHF () real-time time-dependent density functional theories (RT-TDDFT). RG-TDDFT (); octpus-rt-tddft (); gceed (); Meng-Kaxiras-JCP2004-rt-tddft-AO () Here, the electron density operators and associated matrices play central roles in the description of real-time electron dynamics. In particular, compared to the class (a) RT-TDDFT 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 RT-TDDFT. One utilizes time-dependent mutually independent orthogonal orbitals obeying the time-dependent Kohn-Sham (TDKS) equation, Rubio-JCP2004-propagation-tdks () where a time-dependent 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 one-electron 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


where is an orbital label and runs over all occupied orbitals. The Kohn-Sham Fock operator, , depends on the occupied independent Kohn-Sham 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 RT-TDDFT is based on the direct time propagation of the electron density matrix without any propagation of the KS orbitals. Voorithr-rt-TDDFT-density-mat () The KS orbitals are determined through a self-consistent 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 time-dependent Fock matrix, , and light-electron coupling Hamiltonian, ,


Here, . Note that the Fock matrix depends on the time-dependent density matrix at each electron propagation time-step. 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 non-linear 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 time-in-local one and orbital occupation manner for the time dependent KS orbitals as referential independent orbitals. book-fundamental-TDDFT-Gross () 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 donor-acceptor 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 Born-Markov 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 time-dependent manner.

The second approaches are the study of electron transfer dynamics in a test molecular donor-acceptor 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 (multi-layer MCTDH) by Thoss Thoss-TiO2 () and Xie et al. Lan (), respectively. The former author, Thoss, is one of the founders of the multi-layer MCTDH ML-MCTDH () 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 photo-induced electron-transfer process in the alizarin-TiO2 system as a dye-semiconductor 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 four-state electronic Hamiltonian using a GDF representation and DFT calculations and then investigated the multi-electronic-state 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. Scholes-rev-excitons ()

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 (); Thoss-TiO2 () However, here, we have added a density matrix time propagation scheme and light-electron 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 well-known Löwdin orthogonalized atomic basis function, Szabo ()


where the orthogonalized Löwdin atomic orbitals(AO) are expressed by


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 i-th Löwdin atomic orbital, , obtained by Löwdin orthogonalization is the same as that of the original i-th 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


where we used the property of as a real-symmetric matrix. Here, denotes the i-th sub group while is the number of sub groups in the system.

The diagonalization of diagonal block corresponding to subgroup ,


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 spin-orbit couplings. However, because we do not treat spin-orbit 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 ,


which yields the GDF matrix,


In this form, the sub-matrices in the diagonal blocks, , are the diagonal matrices corresponding to the local group eigen energies, while the sub-matrices placed at off-diagonal 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:




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


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 light-electron couplings are , , related to Eq. (9). Here, boldface denotes a vector in three-dimensional Cartesian space, and denotes a composite variable of the electron position in three-dimensional space. The first and second operators are responsible for the light-electron couplings in length and velocity forms. Faisal-LG-VG (); ACP-laser () In this article, we neglect the non-adiabatic 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 Liouville-von Neumann equation and the associated density matrix as follows: where . Here, is a light-electron coupling matrix. The light-electron couplings are described by for the length gauge and for the velocity gauge. Faisal-LG-VG (); ACP-laser () Here, and are the three-dimensional electric field vector and electromagnetic field vector potential, which generally depend on a point in a three-dimensional 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. Faisal-LG-VG (); ACP-laser () 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 fourth-order 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 spin-restricted 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,


For simplicity of presentations in this article, we employed the spin-restricted 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 spin-unrestricted and charged cases will be reported in our future articles.

By using this scheme, we can analyze the dynamical electron re-distribution in molecular aggregate systems driven by the differences in the electron affinities between local groups and the optical-electron 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 donor-acceptor system; the results are used as the first illustration of our scheme. The second example is the unpaired electron dynamics in an excited 20-mer 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 light-harvesting antenna systems LH-anthena-ring (); 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. atom-ring-Ryd-spin (); Eisfeld-ultracold-coins () 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: NPTL-TCNE dimer

Here, by applying the present theoretical method to a NPTL-TCNE dimer system, we observed the initiation of charge separation triggered by the electron-light 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 (); adia-elec-affinity-TCNE (), as is the adiabatic ionization potential of NPTL (8.14 eV) NIST (); ionization-energy-NPTL ()), the NPTL-TCNE set makes an ideal donor-acceptor molecular pair for checking the description of charge transfer with using the present GDF scheme.

Figure 1: (a) Schematics of the relative orientation of the monomers in a NPTL-TCNE dimer. The Cartesian axes are also shown as an eye guide. Atom number labels are also shown. The molecular planes of these two monomers are parallel to the X-Y plane. The central C(23)-C(24) bond line of TCNE is parallel to the Y axis and this is placed on the C(3)-C(6) of the right aromatic ring in NPTL. Note that C(3)-C(6) line is also parallel to the Y-axis. We modulate the distance of these two molecular planes as the difference in their Z-coordinates, namely, the distance of these two molecular planes. (b) HOMO of NPTL(donor) and LUMO of TCNE(acceptor). Orbitals were obtained both for two molecules as isolated systems. at the ab initio level of DFT/6-31Gd/CAMB3LYP. NPTL HOMO and TCNE LUMO energies are -0.2622 (Hartree) and -0.1367, respectively.

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/6-31G(d) level of theory using the CAM-B3LYP exchange correlation functional.CAMB3LYP () The principal axis of NPTL is parallel to the X-axis, while that of TCNE is parallel to the Y-axis. The molecular planes of these flat molecules are parallel to the X-Y plane. Although we treat a dimer here, we partly considered the reported crystal data crystal-NTFL-TCNE-1967 () with respect to the relative orientation as explained in the figure caption. Note that the relative orientation employed here gives rise to a non-vanishing 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 6-31G(d) and functional CAM-B3LYP 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.

Figure 2: Charge migration dynamics in a NPTL-TCNE dimer. The panels include the cases with and without continuum light fields of a wave length of 700 nm and an electronic field strength of 0.02 a.u. The distances between the two monomer planes are 3.0, 3.5 and 4.0 Å  in the left, middle and right columns, respectively. The three panels, (a/b/c), in the top row are the cases without any external light fields. The other nine panels, (d–l), concern the cases with a light field. Normalized polarization vectors corresponding to the 2nd, 3rd and 4th rows are, respectively, , , and .
Figure 3: Charge migration dynamics in a NPTL-TCNE dimer. The panels include the cases with and without continuum light fields of a wave length of 700 nm and an electronic field strength of 0.02 a.u. The distances between the two monomer planes are 3.0, 3.5 and 4.0 Å  in the left, middle and right columns, respectively. The twelve panels, (a–l), concern the cases with a light field. Normalized polarization vectors corresponding to the 1st, 2nd, 3rd and 4th rows are, respectively, , , , and .

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 time-dependent 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 X-axis. 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 field-free 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 back-and-forward 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 Y-axis, we found negligible charge separation, similar to the field-free 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 time-dependent 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 Z-components. First, the results in cases including both X- and Y-components 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 Y-component 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 back-donation caused by the addition of the X-component to the polarization vectors. Interestingly, the Y-polarization of light field in this molecular configuration suppressed the charge separation which can be read from the cases with the combination of Y- and Z-components 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 electronic-bath 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 donor-acceptor system, we obtained the following two findings; (I) a time-dependent 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 donor-acceptor 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: 20-mer 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 ()


where with and being the i-th natural orbital and its corresponding natural population, respectively. Here, denotes a one-electron density operator for the whole system. This treatment enables us to obtain an information concerning the polyradical features of the complex system. b12-nad ()

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 20-mer ethylene

Figure 4: Schematics of the relative positioning of a 20-mer circle of ethylene molecules. The ethylene monomers are identified here by the attached numbers. The Cartesian axes are also shown to guide the eyes of readers.

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 C-C 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 CAM-B3LYP/6-31G(d) level of theory. We also employed the same level of ab-initio calculation for the 20-mer ethylene system as that used for ethylene monomer.

iii.2.3 Exciton migration dynamics

Figure 5: Time-dependent behaviors of quantities of unpaired electrons in the 20-mer circle of ethylene molecules. Unpaired density matrix for the whole system, Eq. (13), was applied to Mulliken population analysis for group monomers throughout (a-i). In the upper three panels, (a/b/c), no group diabatization was employed; that is, the whole system is treated as one monomer and the initial state corresponds to a single electronic excitation from the local HOMO to the LUMO orbital. In contrast, panels (d/e/f) show results of the group diabatic localizations and only 1-st monomer is initially excited from local HOMO to LUMO orbital while 1st and 11th monomers are initially excited in the panels (g/h/i) with the group diabatization procedure. The three panels (a/d/g) correspond to the cases without the light field while the other panels concern light-electron couplings. The parameters of the light fields are as follows; the wavelength and strength are 700 nm and 0.02 atomic unit (a.u.), respectively. and normalized polarization vectors are for (b/e/h) and for (c/f/i). The horizontal and left vertical axes are the time in femtoseconds and the quantity of unpaired electrons for monomers, respectively. The numbers shown on the right vertical axis are the monomer labels. See the main text for details of these plots.

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 20-mer 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 (a-c) 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 NPTL-TCNE dimer. The details of parameters are included in the figure captions. The polarization direction vectors of the light field used here have in-plane and out-of-plane 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 X-Z 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 time-dependent 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 X-Z plane. Thus, the patterns of unpaired electron dynamics is in an inverse relationship between these two cases for the X-Z 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 1-st monomer are separately transferred to 11-th monomer at 10 fs, accompanied by a moderate broadening in the distribution over the monomers. Symmetry with respect to the X-Z 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 X-Z 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 X-Z-reflection 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 two-site single excitations from the local HOMO to the LUMO for the 1-st and 11-th 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 X-Z 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 X-Z plane for these two cases. The difference in the Fourier spectrum and time-dependent 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 two-Frenkel-exciton 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.

Figure 6: Time-dependent Mulliken charge in the 20-mer circle of ethylene molecules. Panels (a/b/c) show results of the group diabatic localizations and only 1-st monomer is initially excited from local HOMO to LUMO orbital while 1st and 11th monomers are initially excited in the panels (d/e/f) with the group diabatization procedure. The three panels (a/d) correspond to the cases without the light field while the other panels concern light-electron couplings. The parameters of the light fields are as follows; the wavelength and, strength are 700 nm, 0.02 atomic unit (a.u.), respectively. and normalized polarization vectors are for (b/e) and for (c/f). The horizontal and left vertical axes are the time in femtoseconds and the scaled unpaired electron quantity for the monomers, respectively. The numbers shown on the right vertical axis are the monomer labels.

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 20-th 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.

  • Compared to the panels of (e), (f), (h) and (i), in Fig. 5 with respect to the light field case with initial excitations, the corresponding panels (b), (c), (e) and (f) of in Fig. 6 show the same shapes of tracks and rectification tread for 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 one-dimensional 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 laser-electron 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 dual-layer 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 (), Pipek-MezeyPM () 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. env-assist-qw-eet () Furthermore, in systems containing metal atoms, the spin-orbit couplings become important for excited state dynamics involved with inter-system crossings. Daniel-ACR-SOC-QM-DYN (); Fedrov-AIMS-SOC () 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 NPTL-TCNE dimer and 20-mer 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.

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, “Next-Generation Supercomputer project” (the K computer project) and “Priority Issue on Post-K computer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use). The part of computations in the present study was performed using Research Center for Computational Science, Okazaki, Japan.


  • (1) V. May and O. Khn, Charge and Energy Transfer Dynamics in Molecular Systems, (Wiley-VCH, 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. Aspuru-Guzik, J. Chem. Phys. 129, 174106 (2008); P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. Aspuru-Guzik, 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 Born-Oppenheimer 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, (Wiley-VCH, 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 Time-Dependent 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 NPTL-TCNE dimer

In Fig. 7, we summarized the mode analysis of GDF electron dynamics in the dimer NPTL-TCNE 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 non-resonant to the typical local electronic transitions because the energy gaps of local HOMO-LUMO 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 y-z 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 HOMO-LUMO 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 20-mer circle of ethylene molecule

Fig. 8 gives the electronic mode analysis in the cases of a 20-mer 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 induced-dipole 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 non-vanishing 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 HOMO-LUMO 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 non-resonant 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 light-field 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 20-mer circle of ethylene molecule

Here we provide the schematic movies of dynamics of unpaired electron in the system of 20-mer ethylene system, as mov-20-mer-two-excitons-1.gif and mov-20-mer-two-excitons-2.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 x-y plane on which all the center of masses (COMs) of ethylene monomers are placed,




the positions of COM of i-th monomer,


and being the quantity of unpaired electrons in i-th 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 20-mer circle of ethylene molecules respectively in Tab. 1, 2 and 3. The associated local orbital energies are also given from the local HOMO-2 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.

Figure 7: Electronic mode analysis in the cases of a dimer NPTL-TNCE system with group diabatic ground state, which are included in the top four raws. The distance of two molecular plane is 4.0 Å. The system is shined by the continuum light wave with 700nm wavelength (corresponding to one photon energy, 0.065 a.u.) and normalized polarization vector of (x,y,z) = , , and , for the panels of (a/f), (b/g), (c/h) and (d/i), respectively. The panels (a/b/c/d) are the Fourier transformations of the time sequence of induced dipole moment of the whole system which are respectively shown in the panels (f/g/h/i). (e) is the Fourier transformation of the continuum light field used in the panels of (d) and (i) while (j) is the electric field in a time domain corresponding to (e).
Figure 8: Mode analysis of electron in the cases of a 20-mer circle of ethylene molecule with initial excitation. The system is shined by the continuum light having 700nm wavelength (corresponding to one photon energy, 0.065 a.u.) and normalized polarization vector of (x,y,z) = . The panels in left column are the Fourier transformations of those placed at just the right position in the figure. The panels (a/c) are the Fourier transformations of the time sequences of induced dipole moment of the whole system for the cases with the initial local excitation of 1st ethylene monomer while (b/d) is that for the initial local excitations of 1st and 11th monomers. The cases of the panels (a/b) associated correspondingly with (f/g) does not include light field while the panels (c/d) corresponding to (h/i) includes. (e) is the Fourier transformation of the continuum light field applied in the cases of the panels (c/d) and correspondingly (h/i). The panels (f/g/h/i) show the induced dipole moment of the whole system corresponding to (a/b/c/d) in this order. The panel (j) displays the electric field of the light in the time domain corresponding to the panel (e).
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
Table 1: Transition dipole moment vector matrix in the group diabatic representation for the NPTL with the distance of molecular planes being 4.0 Å. Unit is atomic unit. Data are shown from the local HOMO-2 ( 32-th local orbital state ) to local LUMO+2 ( 37-th local orbital state ) of this monomer. Note that the origin employed for the evaluation of the electron position is placed at the center of mass of the total dimer system. In this order, local MO energies are Hartree., , , , and . Here, local HOMO and LUMO are labeled by 34 and 35, respectively. As a supporting information, HOMO and LUMO energies in the totally isolated system are -0.2622 Hartree and 0.0115, respectively.
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
Table 2: Transition dipole moment vector matrix in the group diabatic representation for the TCNE with the distance of molecular planes being 4.0 Å. Unit is atomic unit. Data are shown from the local HOMO-2 ( 30-th local orbital state ) to local LUMO+2 ( 35-th local orbital state ) of this monomer. Note that the origin employed for the evaluation of the electron position is placed at the center of mass of the total dimer system. In this order, local MO energies are Hartree, , , , and . Here, local HOMO and LUMO are labeled by 32 and 33, respectively. As a supporting information, HOMO and LUMO energies in the totally isolated system are -0.3901 Hartree and -0.1367,respectively.
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
Table 3: Transition dipole moment vector matrix in the group diabatic representation for the 1st ethylene monomer in the 20-mer ethylene aggregate system. Unit is atomic unit. Data are shown from the local HOMO-2 ( 6-th local orbital state ) to local LUMO+2 ( 11-th local orbital state ) of this monomer. Note that the origin employed for the evaluation of the electron position is placed at the center of mass of the total dimer system. In this order, local MO energies are Hartree, , , , and . Here, local HOMO and LUMO are labeled by 8 and 9, respectively. As a supporting information, HOMO and LUMO energies in the totally isolated system are -0.3277 Hartree and 0.0756,respectively.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description