GW/BSE Nonadiabatic Dynamics Simulations on Excited-State Relaxation Processes of Zinc Phthalocyanine-Fullerene Dyads: Roles of Bridging Chemical Bonds†

In this work, we employ electronic structure calculations and nonadiabatic dynamics simulations based on many-body Green function and BetheSalpeter equation (GW/BSE) methods to study excited-state properties of a zinc phthalocyanine-fullerene (ZnPcC60) dyad with 6-6 and 5-6 configurations. In the former, the initially populated locally excited (LE) state of ZnPc is the lowest S1 state and thus, its subsequent charge separation is relatively slow. In contrast, in the latter, the S1 state is the LE state of C60 while the LE state of ZnPc is much higher in energy. There also exist several charge-transfer (CT) states between the LE states of ZnPc and C60. Thus, one can see apparent charge separation dynamics during excited-state relaxation dynamics from the LE state of ZnPc to that of C60. These points are verified in dynamics simulations. In the first 200 fs, there is a rapid excitation energy transfer from ZnPc to C60, followed by an ultrafast charge separation to form a CT intermediate state. This process is mainly driven by hole transfer from C60 to ZnPc. The present work demonstrates that different bonding patterns (i.e. 5-6 and 6-6) of the C−N linker can be used to tune excited-state properties and thereto optoelectronic properties of covalently bonded ZnPc-C60 dyads. Methodologically, it is proven that combined GW/BSE nonadiabatic dynamics method is a practical and reliable tool for exploring photoinduced dynamics of nonperiodic dyads, organometallic molecules, quantum dots, nanoclusters, etc.

In this work, we employ electronic structure calculations and nonadiabatic dynamics simulations based on many-body Green function and Bethe-Salpeter equation (GW/BSE) methods to study excited-state properties of a zinc phthalocyanine-fullerene (ZnPc-C 60 ) dyad with 6-6 and 5-6 configurations. In the former, the initially populated locally excited (LE) state of ZnPc is the lowest S 1 state and thus, its subsequent charge separation is relatively slow. In contrast, in the latter, the S 1 state is the LE state of C 60 while the LE state of ZnPc is much higher in energy. There also exist several charge-transfer (CT) states between the LE states of ZnPc and C 60 . Thus, one can see apparent charge separation dynamics during excited-state relaxation dynamics from the LE state of ZnPc to that of C 60 . These points are verified in dynamics simulations. In the first 200 fs, there is a rapid excitation energy transfer from ZnPc to C 60 , followed by an ultrafast charge separation to form a CT intermediate state. This process is mainly driven by hole transfer from C 60 to ZnPc. The present work demonstrates that different bonding patterns (i.e. 5-6 and 6-6) of the C−N linker can be used to tune excited-state properties and thereto optoelectronic properties of covalently bonded ZnPc-C 60 dyads. Methodologically, it is proven that combined GW/BSE nonadiabatic dynamics method is a practical and reliable tool for exploring photoinduced dynamics of nonperiodic dyads, organometallic molecules, quantum dots, nanoclusters, etc.

I. INTRODUCTION
Donor-acceptor dyads have attracted a broad range of research interests owing to their successful applications in photo-fuel production, photo-electricity conversion, optoelectronic devices, etc. [1][2][3][4]. Their optoelectric performances heavily depend on initial photoinduced charge separation and subsequent carrier dynamics. Clearly, correctly understanding these dynami-cal processes is very important for rationally regulating and designing new types of donor-acceptor dyads with superior optoelectrical properties. Thus, numerous experimental and theoretical works have been carried out in the past decades to study excited-state properties of dyads with distinct kinds of donor-acceptor interfaces [5][6][7][8][9][10][11].
One kind of the most popular donor-acceptor dyads is formed by combining both zinc phthalocyanine (ZnPc) molecules and fullerenes (C 60 ), e.g. ZnPc-C 60 , etc. These ZnPc-C 60 complexes are overall categorized into two types according to whether they are formed by chemical bonds or not. Non-covalently bonded ZnPc-C 60 complexes are primarily formed through van der Waals interaction. Alternatively, these ZnPc-C 60 dyads can be generated via chemical bonds connecting both donor and acceptor fragments. In these complexes, ZnPc and C 60 usually act as electron donor and acceptor, respectively, and their unique optoelectronic properties have been extensively studied in the past years [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Once exciting the donor fragment, i.e. ZnPc, singlet excitons are first generated and then migrate toward ZnPc-C 60 interfaces near which a series of nonadiabatic transitions such as internal conversion and intersystem crossing processes could occur. Finally, these excitons are converted into charge-transfer excitons followed by charge separation to generate free electrons and holes. Furthermore, these photoinduced interfacial dynamical processes are very complicated and seriously dependent on relative spatial orientations of ZnPc and C 60 in noncovalently bonded dyads, and chemicalbond properties connecting ZnPc and C 60 in covalently bonded dyads, as observed in many experimental and theoretical studies [22, 24-26, 28, 29]. In order to uncover chemical bonding effects on excited-state properties and charge separation dynamics of covalently bonded ZnPc-C 60 dyads, several ZnPc-C 60 dyads have been synthesized and their photochemical and photophysical properties have been explored [30][31][32][33][34]. Nevertheless, mechanistic details on excited-state dynamics, in particular effects of chemical bonds connecting both donor and acceptor fragments, remain elusive.
Theoretically, there are few works carried out to address excited-state properties and photo-driven dynamical processes of ZnPc-C 60 dyads. Santos and Wang carried out Ehrenfest dynamics simulations to explore photoinduced electron and hole transfer dynamics between ZnPc and C 60 in two different covalently bonded ZnPc-C 60 complexes [29]. Their results demonstrated that chemical bonds connecting ZnPc and C 60 have remarkable influences on photoinduced charge separation dynamics. However, molecular-orbital-based mean-field methods are not suitable for simulating excitation energy transfer due to lacking of electron-hole interactions. Later, our theoretical works based on the linearresponse time-dependent density functional theory (LR-TDDFT) method revealed that excited-state electronic structures and relaxation dynamics of ZnPc-C 60 dyads are significantly regulated by spatial orientation in noncovalently bonded dyads and chemical bonding patterns in covalently bonded ones [28,35]. However, in previous works, only C−O linking groups are considered. Whether other chemical linkers, such as C−N groups, have similar effects, is not clear. Motivated by this question, we here developed and employed GW/BSE (combined many-body Green function and Bethe-Salpeter equation) nonadiabatic dynamics simulations to study excited-state properties and photoinduced dynamics of ZnPc-C 60 dyads linked by two C−N bonds but with different chemical bonding patterns.

A. GW and BSE methods
The GW and BSE methods are derived from manybody perturbation theory (MBPT) [36][37][38][39][40][41][42]. The core object in MBPT is time-ordered Green function G, which is expressed as where Φ represents a ground-state wavefunction and T is a time-ordering operator.ψ † (x ′ ,t ′ ) andψ(x, t) represent electron creation and annihilation operators in the Heisenberg picture. t represents time and x=(r, σ) are space and spin coordinates (r: space; σ: spin). G describes time-dependent propagation of electrons and holes and only depends on x and t. Meantime, , depends on frequency ω. One-particle Green function G σ (r, r ′ , ω) could be expressed with one-electron wavefunctions ψ σ m (r) and energies ε σ m calculated in meanfield potential of electrons where indexes i and a represent occupied and unoccupied one-electron states with spin number σ, and η is a small positive number. Non-interacting polarization χ 0 is related to oneparticle Green function G σ (r, r ′ , ω) and dielectric function ε(r, r ′ , ω) with random-phase approximation (RPA) is finally given by in which v represents Coulomb interaction. With ε, one can get screened Coulomb interaction W (r, r ′ , ω) as then, self-energy is written as where e iηω1 is merely used to enforce a correct timeorder form of self-energy. With these definitions, GW quasiparticle energies ε GWσ n can be obtain by which is also called one-shot G 0 W 0 scheme. In addition, one can also further iterate GW equations by updating Eqs.(2−7) with new quasiparticle eigenvalues, while electron wavefunctions are kept fixed without update in this procedure. This procedure is called eigenvalue-only self-consistent GW, i.e. evGW or G n W n in literature, in which n means iteration cycles for Green function and screened Coulomb interaction. Of course, there are GW methods with updated eigenvalues and wavefunctions, however, which are beyond the scope of our present work. Interested readers are referred to recent literature [43][44][45][46][47][48][49].
On the basis of GW results, BSE can be used to calculate excitation energies, which is expressed as where L 0 represents the non-interacting counterpart of L and the electron-hole interaction kernel K(r 3 , r 5 ; r 4 , r 6 ) is written as Then, after further mathematical derivation, one can get excitation and de-excitation amplitudes of electronhole transition pairs for the lth excited state which is similar to that of LR-TDDFT method, i.e. the Casida equation [50]. Thus, BSE excitation energies can be similarly obtained by in which Ω l represents excitation energy of the lth excited state, and X l and Y l are the corresponding excitation and de-excitation amplitudes. Matrix elements of A and B are expressed as (13) and in which i, j represent occupied states while a, b represent unoccupied states. α S/T is a constant which equals to 0 for triplet states and 2 for singlet states.

B. Fewest-switches surface-hopping method
Nonadiabatic molecular dynamics simulations play important roles in understanding ultrafast photochemical and photophysical processes of molecules, biological systems, and materials. Nowadays, there are a lot of simulation methods proposed and widely applied [10,11,. One kind of the most popular simulation algorithms is trajectory-based surface hopping methods [73,74]. The main idea of surface hopping method is that nuclei propagate on a specific adiabatic potential energy surface and have probabilities of hopping to other adiabatic ones. The fewest-switches surface hopping method, proposed by Tully et al. [73,74] is one of the most popular surface hopping methods. Recently, we have developed and implemented this surface hopping simulation method at the LR-TDDFT level, and successfully applied this method to simulate many complex photoinduced processes [28,35,[75][76][77][78]. However, sometimes, TDDFT calculated excited-state electronic properties are heavily dependent on the used density functionals. They could vary qualitatively from one to another. By contrast, the GW/BSE method based on MBPT is more robust, in particular evGW/BSE, which provides reliable results of excited states. Thus, in this work, we develop and implement the fewest-switches surface hopping method combined with the GW/BSE method.
Here we give a brief introduction. The timedependent Schrödinger equation is written as whereĤ 0 (R, r) is zero-order electronic Hamiltonian while r and R represent spatial coordinates of electron and nuclear, respectively. Here, time-dependent total electronic wavefunction ψ(R, r, t) is expanded as where c i (t) is a time-dependent expansion coefficient of the ith electronic state and Φ i (R, r) is corresponding adiabatic electronic wavefunction at nuclear configuration R. Taking Eq. (17) into Eq.(16) and multiplying ⟨Φ i (R, r)| from left side, one can get in which d ij and v(t) are adiabatic derivative couplings and nuclear velocities, respectively. E j is the eigenvalue ofĤ 0 (R, r) with corresponding eigenfunction Φ j (R, r). Eq. (18) is the final electronic propagation equation. Relevant fewest-switches criterion that computes electronic transition probability from i to j states, which is written as

C. Numerical algorithm for nonadiabatic couplings
Numerical algorithm for calculating nonadiabatic couplings has been developed and implemented by us at the LR-TDDFT level [28,35,[75][76][77][78]. Since BSE is similar to LR-TDDFT in mathematical forms, it is easy and straightforward to modify the algorithm to calculate nonadiabatic couplings at the GW/BSE level. Similarly, the Kth excited-state wavefunction Φ K at the BSE level is written as linear combination of GW molecular orbitals ψ a i : where ω stands for a linear combination coefficients of ψ a i , and i, a represent occupied states and unoccupied states, respectively. Nonadiabatic couplings τ KJ between states K and J can be written as DOI:10.1063/1674-0068/cjcp2109162 c ⃝2021 Chinese Physical Society in which i, j represent occupied states while a, b represent unoccupied states. Note that mathematical expression of τ KJ within the GW/BSE framework is the same as that within the LR-TDDFT one [75,79]. After mathematical derivation, the final expression of τ KJ can be written as where P ij is an additional phase factor which depends on the ordering convention of molecular orbitals. Time differentiation on molecular orbitals can be obtained by finite-difference scheme (23) in which ϕ p (t) and ϕ q (t+∆t) represent MOs at time t and t+∆t, respectively. Because our used GW/BSE calculations do not involve updating molecular orbitals, molecular orbitals from converged DFT calculations are used in GW/BSE calculations as ϕ p (t) and ϕ q (t+∆t) in Eq.(23).

D. Simulation details
Geometries of the ZnPc-C 60 dyad were optimized using the B3LYP+D3 method [80][81][82][83][84]. The C, H, O, and N atoms were described with the cc-pVDZ basis sets [85]. For the Zn atom, LANL2DZ basis sets were employed to describe outer-valence electrons while innercore electrons were treated with pseudopotential [86]. Geometry optimization was performed with Gaussian 16 [87]. Energy decomposition analyses (EDA) were performed at the B3LYP+D3/TZP level of theory with ADF2018 [88,89]. Ground-state molecular dynamics simulations were performed at the PBE level with DZVP-MOLOPT-SR-GTH basis sets and Goedecker-Teter-Hutter pseudopotentials [90][91][92][93][94], which were implemented in the QUICKSTEP module of CP2K-5.1 [95][96][97]. Starting from optimized structures, we first heated the system to 300 K and equilibrated it for about 1 ps with a 1 fs time step, during which the Nosé-Hoover chain thermostat technique (chain length: five) was employed [98,99]. After that, 2 ps microcanonical dynamics simulations were performed. All GW and BSE calculations are performed with MOLGW-2.C [100] with cc-pVDZ basis sets used for all atoms except cc-pVDZ-PP basis sets and pseudopotentials for the Zn atom [85,101]. In evGW/BSE calculations, GW calculations with two iteration steps were found to give converged results in Table S2 in Supplementary materials, which was thus employed in all the evGW/BSE calculations. The required molecular orbitals and initial eigenvalues were calculated with the DFT method. All nonadiabatic dynamics simulations were conducted using our developed GTSH package [59]. Electronic transition density analyses on evGW/BSE results were calculated using MULTIWFN3.6 [102,103].

III. RESULTS AND DISCUSSION
As the first step to decipher excited-state properties of the C-N linked ZnPc-C 60 complex, we first optimize the two ground-state structures with different linking patterns, i.e. 5-6 and 6-6 types, at the B3LYP+D3 level. In the former, two N atoms of ZnPc are bonded to two C atoms shared by one pentagonal and one hexagonal rings of C 60 ; while, in the latter, two N atoms are bonded to two C atoms shared by two hexagonal rings of C 60 (see FIG. 1). As can be seen, these two structures are overall similar to the previously studied C−O linked 5-6 and 6-6 configurations. Additionally, the energy decomposition analysis based on natural orbitals from chemical valence method is performed with respect to both ZnPc and C 60 fragments. The results are listed in Table I. Total interaction energy E tot is composed of four kinds of contributions, i.e. E tot =E Pauli +E ele +E orb +E dis , in which E Pauli is exchange repulsion energy due to Pauli principle among fragments, E ele represents electrostatic interaction energy among fragmental charge densities, E orb is energy due to fragmental orbital mixing, and E dis is dispersion correction energy among fragments (see more details and discussion in Supplementary materials). As with the C−O linked ZnPc-C 60 dyad [35], E tot of the 6-6 configuration is larger than that of the 5-6 configuration (116.1 vs. 93.6 kcal/mol), which indicates that the 6-6 configuration is more stable than the 5-6 one. This difference is mainly originated from E orb , which is −540.2 kcal/mol for the 6-6 configuration and −517.3 kcal/mol for the 5-6 configuration. In comparison, the other three kinds of interactions (i.e. E Pauli , E ele , and E dis ), contribute less to the difference.
To elucidate excited-state properties of these two ZnPc-C 60 complexes, we perform evGW/BSE calcu-  lations to obtain their excitation energies and electronic transition densities of the lowest 10 singlet excited states, as shown in FIG. 2, FIG. S2 and Table S3 in Supplementary materials. As discussed above, molecular orbitals from converged DFT calculations are input as initial wavefunctions and will not change in evGW/BSE calculations. In order to check the influence of exchange-correlation functionals on evGW/BSE results, we have tested several exchangecorrelation functionals and found that evGW/BSE calculated excited-state properties are not sensitive to the functionals used (see Tables S4 and S5 in Supplementary materials). Therefore, the B3LYP functional is chosen in all the evGW/BSE calculations. In Table S3, one can see that excitation energies from S 2 to S 10 for the 5-6 configuration range from 1.649 eV to 2.131 eV and energy difference between all the nearby states is less than 0.15 eV (i.e. 3.5 kcal/mol). Among these states, the S 2 state is of pure CT character and its electron is localized on C 60 while hole is localized on ZnPc. The S 7 and S 8 states are of mixed LE and CT character (see Table S6 in Supplementary materials). All the other states are pure LE states of either C 60 or ZnPc (see Table S3 in Supplementary materials). In addition, oscillator strengths of the lowest six excited states are very small. In contrast, excited states from S 7 to S 10 have larger oscillator strengths. As mentioned before, ZnPc is usually chosen as electron donor and the lowest LE state on ZnPc is chosen as initial populated singlet state in the present non-adiabatic dynamics simulations.
As to the 6-6 bonding configuration, the LE state of ZnPc is the lowest singlet excited state indicating that both electron and hole localize on ZnPc (see Table S3 in supplementary materials). Therefore, if this state is populated in the Franck-Condon region, there is small probability for ultrafast charge separation taking place. The same situation is also seen in the C−O bonded ZnPc-C 60 dyad studied recently [35]. According to these results, we will next focus on photoinduced dynamics of the 5-6 bonding configuration.
Despite valuable results from electronic structure calculations, photoinduced dynamical processes of the present ZnPc-C 60 dyad remain unclear. In order to figure out the details of photodynamics, we carried out nonadiabatic dynamics simulations based on the fewestswitches surface hopping method at the evGW/BSE level. The classical path approximation is widely used for simulating photoinduced processes without large conformational changes or bond closing/breaking [28,35,[59][60][61][62][75][76][77][78]. In dynamics simulations, the lowest bright LE state of ZnPc is chosen as the initial state, 300 initial conditions are randomly chosen from a pre-defined microcanonical trajectory; for each initial condition, 200 surface hopping trajectories are propagated for 500 fs. The results are finally averaged over 300×200 trajectories.
We first plot time-dependent populations of the lowest 10 excited states, as shown in panel of FIG. 3(a), in which weights from S 4 to S 10 are summed together since the LE state of ZnPc varies among these states due to their close energies at different initial conditions. It is clear that these excited states population decays to zero in an ultrafast way within ca. 400 fs. Fitting the population curve with a simple mono-exponential function gives an effective excited-state population decay rate constant of 159 fs (S 4 to S 10 ). At the same time, the S 3 weight first increases to a maximum of 0.35 at 150 fs and then decays to nearly zero at 500 fs. The S 3 population rising and decay rate constants are predicted to be 116 and 193 fs, respectively. In contrast, both S 1 and S 2 populations increase monotonously within 500 fs and do not reach their maxima. The S 2 weight increases to 0.85 while that of the S 1 state increases slowly to 0.15 at 500 fs (see FIG. 3(a)). In a similar way, their population rising rate constants are estimated to be 223 and 1386 fs, respectively. Clear analysis reveals that the ultrafast decay of the S i (i=4−10) state to the S 3 state is due to their extremely small energy differences and large nonadiabatic coupling terms, as shown in FIG. 3(c) and FIG. 4, and Table S7 in supplementary materials. As can be seen, most energy gaps between E i+1 −E i (i≥3) are less than 5 kcal/mol, and the average energy gap is 1 kcal/mol (see FIG. 3(c)). On the other hand, nonadiabatic coupling terms between i+1 and i states (i≥3) are very strong and average values are larger than 75 ps −1 (see FIG. 4). Once arriving at the S 3 state, a little larger energy gap between S 3 and S 2 , an average value of 2.5 kcal/mol, delays further nonadiabatic transition to the S 2 state (see FIG. 3(b)). Nevertheless, it still results in a fast decay of the S 3 state in 400 fs because of still small energy gap. Finally, unlike the fast population of the S 2 state, its de-population dynamics to the S 1 state is very slow due to large energy gaps between them and Since different CT and LE (ZnPc or C 60 ) electronic states are involved in nonadiabatic dynamics simulations, nonadiabatic transitions will drive exciton, electron, and hole transfer among ZnPc and C 60 fragments. With the fragment-based exciton analysis method proposed in our previous works [28,35,75,76,78], we have plotted time-dependent electron and hole amounts on both ZnPc and C 60 fragments in FIG. 5. One can see that at the beginning of nonadiabatic dynamics simulations, electron and hole are mainly located on ZnPc with weights of 0.8, which agrees with the |ZnPc * -C 60 ⟩ exciton weight at the initial simulation time (see FIG. 6(a)). Then, both electron and hole amounts located on ZnPc sharply decreases and becomes stable at 200 fs. The corresponding electron and hole population decay rate constants are estimated to be 60 and 64 fs, respectively, according to fitting the curve with two exponential functions. In this process, electron and hole transfer are completely synchronous and in the same direction from ZnPc to C 60 , which is consistent with the feature of excitation energy transfer. This is also seconded by FIG.  6 (a) and (b) where the |ZnPc * -C 60 ⟩ exciton decrease is accompanied with the |ZnPc-C 60 * ⟩ exciton increase (ca. 100 fs). The corresponding population decay and rising rate constants are fitted to be 60 and 54 fs, which agree very well with the above hole and electron population decay rate constant of 60 fs. Afterwards, the remaining electron amount on ZnPc is gradually transferred to C 60 until the end of simulations (see FIG. 5(a) and (b)); while, the hole amount, already transferred from ZnPc to C 60 , starts to move back to ZnPc at 100 fs, because of exciton transfer from |ZnPc-C 60 * ⟩ to |ZnPc + -C 60 − ⟩.
Its hole population then becomes stable at 400 fs with a rising rate constant of 111 fs. At 500 fs, both ZnPc and C 60 possess equal hole amounts revealing a mixture of |ZnPc-C 60 * ⟩ and |ZnPc + -C 60 − ⟩ (see FIG. 6(b) and (c)).
Then, one must be noted that due to large energy gaps between S 2 and S 1 , there are comparable S 2 and S 1 populations at the end of the 500 fs simulations (see FIG. 3(a)). Nevertheless, according to electronic structure calculations, further nonadiabatic decay to S 1 is thermodynamically allowed and must happen but with a longer time. Therefore, one can safely postulate that the hole localized on ZnPc will finally transfer back to C 60 and the pure LE exciton |ZnPc-C 60 * ⟩ will be eventually populated. Interestingly, in the entire dynamics simulations, there is no |ZnPc − -C 60 + ⟩ charge-transfer exciton from C 60 to ZnPc observed, which is consistent with the preceding electronic structure calculations. As shown in FIG. 2, there is no such kind of CT states from C 60 to ZnPc available in the low-lying excited states.
Finally, time-dependent exciton size is useful for a better understanding of time-dependent excited properties of dyads. At the beginning of nonadiabatic dynamics simulations, exciton size is 5.84Å. It is nearly not changed, 5.60Å in the first energy transfer process of 50 fs, due to the local character of both |ZnPc * -C 60 ⟩ and |ZnPc-C 60 * ⟩ excitons. After that, the |ZnPc-C 60 * ⟩ exciton with the LE character hops to the |ZnPc + -C 60 − ⟩ exciton with the CT character, which results in a significant increase of exciton size from 5.60Å to 8.58Å at 500 fs due to charge separation.
The excited-state dynamics of the C−N bonded ZnPc-C 60 dyad is generally similar to that of the C−O bonded ZnPc-C 60 dyad. For example, for the 6-6 bonding configuration, the LE state of ZnPc is the lowest S 1 state [35]. Thus, initial photoexcitation populates this S 1 state and its further charge separation is relatively slow and cannot complete within several hundreds of femtosecond, as demonstrated by corresponding dynamics simulations. However, for the 5-6 bonding configuration, the S 1 state is the LE state of C 60 with very small oscillator strength; while, that of ZnPc is much higher in energy. Importantly, there are some CT states between the LE states of ZnPc and C 60 . As a result, one can see obvious charge separation process during excited-state relaxation processes from higher  [35]. In addition, the energy gap between S 2 and S 1 is relatively smaller in the C−O bonded complex than that in the C−N bonded one. Thus, at 500 fs, the S 1 population is close to 0.7 in the former; while, it is only 0.15 in the latter. This difference makes the exciton size first increase and then decrease in the C−O bonded complex. But, the decrease process is not seen within 500 fs in the C−N bonded complex (see FIG. 7). In one word, the different properties of chemical bonds connecting both ZnPc and C 60 of the ZnPc-C 60 dyad have remarkable influences on its excited-state properties and nonadiabatic dynamics, and thereto optoelectronic properties.

IV. CONCLUSION
In this work we have studied excited-state properties and photoinduced dynamics of both 5-6 and 6-6 bonding configurations of a ZnPc-C 60 dyad with accurate evGW/BSE calculations and related nonadiabatic dynamics simulations. It is found that different chemical bonding patterns between ZnPc and C 60 have remarkable influences on excited-state properties and subsequent nonadiabatic dynamics. In the 6-6 bonding configuration, the S 1 state is the LE state of ZnPc with large oscillator strength and thus is initially populated excited state. However, it is difficult for further charge separation from this S 1 state. In stark contrast, in the 5-6 bonding configuration, the LE state of ZnPc is much higher, which is the S 9 state at the Franck-Condon point. Below this bright state, there are several excited states with obvious CT character between ZnPc and C 60 , and LE character on either ZnPc or C 60 . Therefore, thermodynamically, it is allowed from the higher S 9 to S 1 states via several intermediate states. These points are supported by 500 fs dynamics simulations. In the first 50 fs, one can see an ultrafast excitation energy transfer from ZnPc to C 60 (concerted electron-hole pair transfer). This process is followed by a fast charge separation driven by hole migration back to ZnPc from C 60 within 500 fs. The last nonadiabatic transition to the S 1 state is not observed due to limited simulation time, but is allowed thermodynamically in terms of electronic structure calculations. In combination our previous work [35], we have proven that different chemical bonding patterns and chemical linkers can be used to tune excited-state properties and optoelectronic prop-erties of covalently bonded ZnPc-C 60 dyads. Therefore, if new ZnPc-C 60 like dyads with ultrafast energy transfer properties are expected, the 5-6 bonding configuration might be more preferred. Methodologically, we demonstrate that the combined evGW/BSE nonadiabatic surface-hopping dynamics method is reliable and accurate, which will encourage more theoretical simulations on ultrafast photophysical and photochemical processes of molecules, clusters, etc.
Supplementary materials: Additional theoretical methods, evGW/BSE results with different DFT wavefunctions, evGW/BSE results with different iterative cycles, additional data for averaged NAC and for the 6-6 bonding configuration, and Cartesian coordinates are available.