Computing intramolecular charge and energy transfer rates using optimal modes

One of the most important and fundamental processes in chemical dynamics is that of energy and charge transfer between molecular species. The mechanism itself is highly quantum mechanical in nature and involves strong coupling between nuclear and electronic degrees of freedom. The seminal model for the calculation of the transfer rate for this process was developed by Marcus in the 1950s,^1–3 

which relates the Fermi golden-rule transition rate, k_Marcus to the thermodynamic driving force ΔG^o and the reorganization energy λ. One of the most striking predictions of this theory is that as the free energy difference between the final and initial state increases, the transition rate reaches its maximum, and further increases in the energy off-set leads to slower transition rates. This “inverted regime” was confirmed by Closs and Miller^4–6 nearly 30 yr later.

Numerous improvements to the Marcus model have been presented over the years and a detailed exploration of each is well beyond the scope of this work. Aside from the electronic coupling, V_ab, the other terms in the expression of rate constant arise from a semi-classical approximation to the overlap between the initial vibrational states of the donor and the final states of the acceptor state. It suffices to say that most of the efforts focus upon incorporating more molecular-level details into the couplings and vibrational modes used in computing these terms.

A number of years ago, we developed a time-convolutionless master equation (TCLME) approach for computing state-to-state rate in which the coupling between states depends upon the nuclear coordinates.⁷ This approach incorporates a fully quantum mechanical treatment of both the nuclear and electronic degrees of freedom and recovers the Marcus rate equation in the semiclassical limit. The model itself is parameterized by the vibrational normal mode frequencies, the electronic energies, and energy derivatives at a reference configuration obtained from ab initio quantum chemistry computations. The approach has been intensively used and testified by our group to compute state-to-state transition rates in semi-empirical models for organic semiconducting light-emitting diode and photovoltaics.^8–11 Central to our work here is the use of a mode-projection scheme which parses out a reduced set of nuclear motions primarily coupled to the quantum transition. We refer to the modes identified early in the iterative process as the “primary Lanczos modes” or PLMs.

Using the TCLME approach, we recently investigated triplet-triplet excitation energy transfer between a naphthalene donor and a benzaldehyde (BZ) acceptor linked by a variety of bridging units as shown in Fig. 1.¹² Such systems were the basis of a classic series of experiments by Closs and Miller^4–6 which verified the existence of the Marcus inverted regime. In Ref. 12, we showed that both the autocorrelation function of the electronic coupling operators and the total transfer rate constant calculated using TCLME approach with only PLMs provide an excellent approximation to the exact correlation functions and rate constants computed using all normal modes, as well as to the experimental and recent theoretical values for the transfer rates.¹³ 

In this paper, we analyze the symmetry of PLMs in model donor-bridge-acceptor systems as projected onto the local vibrational modes of the donor (benzaldehyde) and acceptor (naphthalene) moieties. By constructing the local modes from the dominant normal modes in the projection, one can deduce a relation between the PLMs, fragment modes, and normal modes for the entire molecule. The structure of this paper is as follows. We begin by reviewing the TCLME approach and how one generates the PLMs. We then discuss the symmetry of PLMs and their relation to a linear coordinate linking the initial/donor geometry to the final/acceptor geometry of the molecule.

Our approach for constructing the PLMs was presented in Ref. 12. To be self-contained, we give a brief overview of their construction and connection to electronic structure models.

We start with the generic form for the diabatic Hamiltonian, describing two electronic states coupled with N normal modes,

Here, the first term contains the electronic energies, ϵ₁ and ϵ₂, computed at a reference geometry–typically that of the donor or acceptor state. V_ij is the diabatic coupling between them. The second term represents the linear coupling between the electronic and nuclear degrees of freedom given in terms of the mass-weighted normal coordinates q. The diagonal terms give the diabatic displacement forces between the reference geometry and the two states. The remaining two terms correspond to the harmonic motions of the nuclear normal modes, given here in mass-weighted normal coordinates. In the normal mode basis, Ω is diagonal with elements corresponding to the normal mode frequencies, $ω j 2$⁠.

If we choose one of the states as the reference state, then either g₁₁ or g₂₂ will vanish. We further assume Condon approximation to neglect g₁₂ and g₂₁. Now, the Hamiltonian is simplified to

where H_osc is the harmonic oscillator Hamiltonian for the vibrational normal modes.

To calculate transition rate, we transform to adiabatic representation. It is done via Edmiston-Ruendenberg (ER) localization¹³ implemented in Q-Chem and the electronic Hamiltonians on adiabatic and diabatic bases are related by the mixing angle θ,

The diabatic coupling is then given by

The full adiabatic Hamiltonian reads

where the electronic part is in the adiabatic basis with eigenenergies of E₁ and E₂ at certain reference geometry. The nuclear normal modes are described by H_osc, whereas θ is the mixing angle between two states.

To parametrize our Hamiltonian, we assume the diabatics are a good approximation to the actual adiabatic potentials. This allows us to use the gradients of the adiabatic potentials to approximate the diabatic potentials. This approximation is valid when the adiabatic (and diabatic) energy minima are far enough away from the crossing points and the mixing angles between the diabatic and adiabatic states are small. Thus, all parameters in Eq. (6) can be obtained from standard quantum chemical computations. As in our previous work, we use the Q-Chem 4.0 package to obtain the vertical energies from single point CI(S) calculations with 6-31G** basis set at a given reference geometry. We then project the energy gradients onto the vibrational normal coordinates.

Our group has shown that the dynamics of Hamiltonian in the form of Eq. (6) can be expressed in a time-convolutionless mater equation.⁷ The derivation of the TCLME and its Fermi golden-rule limit is lengthy and we give a brief recapitulation of the derivation and motivation next.

We start from a more general Hamiltonian,

then perform a polaron transform using the unitary transformation,^7,14,15

under which the transformed Hamiltonian is written in terms of the diagonal elements

with the renormalized electronic energies,

and off-diagonal terms,

After lengthy derivation we get the autocorrelation function in the Heisenberg representation

where 〈⋯〉 denotes a thermal average over the vibrational degrees of freedom. Using this approach, the golden-rule rates are given by

In a practical sense, we take τ to be some finite time at which the autocorrelation function C(t) has relaxed to zero.

Later, we will use the C_nm(t) to benchmark the convergence of our model with respect to the number of nuclear modes included in constructing the electronic coupling operator, V_nm(t). For our purposes here, an “exact” calculation involves including all nuclear vibrational modes. In our previous work, we showed that both C(t) and the total transfer rate constant, k_nm, calculated using only the first few projected modes provide an excellent agreement with the exact quantities computed using the full set of normal modes, as well as the experimental rates, when parameterized using accurate quantum chemical data.¹² 

The linear coupling parameter g₂₂ in Eq. (6) defines a force along a vector connecting the initial and final equilibrium geometries of the molecule. This vector, along with the diabatic mixing angle can be obtained from quantum chemistry using the ER localization scheme to approximate the donor and acceptor states. By analyzing this force, we can gain the insight into the dynamics of the transition as well as open an avenue for developing improved approximations for transition rates. In our previous work,¹² we presented a Lanczos-base ranking algorithm that projects out a series of nuclear displacements that are most important for the transition. We refer to the highest-ranked mode identified by the algorithm as the “PLM.” The process of finding the PLMs is initiated by defining the vector v₁ = g₂₂. At each step indexed by k, we define a projection operator

and its complement Q_k = I − P_k. We also construct

as the total projection operator for all k ≤ N modes. We then project the Hessian matrix Ω into each subspace viz.

and diagonalize each to obtain eigenvalues and eigenvectors {α_p, M_p} and {α_q, M_q}, respectively. Ω_p and Ω_q are N × N matrices. The first set of these will have a single non-trivial eigenvalue and the second set will have N − k non-trivial eigenvalues. We then collect the non-trivial eigenvectors associated with each to form the orthogonal transformation matrix

and again transform the full Hessian Ω into this new vector space to form the N × N matrix Ω′. At each step in the iteration, the transformed Hessian, Ω′, is in the form of a k × k tri-diagonal submatrix in the upper-left part of the matrix and a diagonal submatrix in the lower-right.

In Fig. 2, we show the PLMs for the intramolecular triplet energy transfer within two representative donor-bridge-acceptor molecules from our study, termed c14ee and d26ae. In both cases, triplet energy is transferred from a BZ donor group to a naphthyl-acceptor group. Several observations can be made from Fig. 2 concerning the primary Lanczos modes for these donor-bridge-acceptor systems. First, there is negligible contribution from the bridging unit. This is not surprising since the transitions involve electronic energy transfer between π-orbitals localized on the donor and acceptor moieties. Since the bridges in these cases are not conjugated, there is very weak electronic coupling between the D/A groups and the bridge itself. Hence, to a good approximation, the electronic contribution from the bridge can be effectively ignored, it simply serves to hold the donor and acceptor in fixed relative positions. Second, the displacement vectors on donor and acceptor moieties are very much akin to totally symmetric normal modes of the corresponding moieties and do not correspond any single normal vibrational mode of the whole molecule.

In Secs. III and IV, we explore the role of symmetry in detail. We first project the PLMs onto the normal modes of both the donor or acceptor moieties and onto entire molecules (donor-bridge-acceptor). Using these reduced PLMs, we compute rate constants and compare to the numerically exact results obtained using the full modes. Then, we do a similar set of projections using simply the geometric change (GC) from initial to final configurations to explore the relation between the PLMs and geometry.

Fig. 3 shows the projection of the PLMs onto the local vibrational modes of D/A moieties in molecules of c14ee and d26ae. In each case, we label the dominant contributions by the irreducible representation (IR) of the moiety. The normal mode basis is generated as follows: First, we take the geometry of donor or acceptor moiety and add the necessary hydrogens. Then, the D_2h symmetry is applied to naphthalene using GaussView 5. The benzaldehyde moiety is unchanged because it has very low symmetry: C_s if planar and C₁ other wise. If one treats the aldehyde group as a local site, one can classify the ring-motions of benzaldehyde in terms of the C_2v irreducible representations. As a result, the assigned irreducible representations are exact for naphthalene but approximate for benzaldehyde. We then assign a weight to each mode viz.

where $M i ⃗$ is the displacement vector of the ith normal mode. The symmetrization of naphthalene introduces another problem in that the projection depends on the orientation of the basis vectors. Thus, we orient the original and symmetrized frames by first aligning the center of masses and then minimize the RMS displacement between atoms. Analyzing the projections, it becomes clear that in all cases, the dominant components belong to totally symmetric A_g or A₁ irreducible representations, in agreement with the intuitive observations from Fig. 2.

This behavior is not unique for the napthyl-bridge-benzaldehyde system. To see this, we repeated the analysis using various acceptor groups as shown in Fig. 4, together with their PLMs and correlation functions for different numbers of projected modes. For both molecules, with 4 projected modes or more, the excellent agreement with exact (all modes) result out to 10 fs indicates that the projection scheme provides an accurate description of the electronic coupling in these cases as well. Moreover, the PLMs are qualitatively similar to the PLM of unchanged c14ee. The displacements on the benzaldehyde moiety are almost identical, while on the phenyl and anthracene moieties, they belong to totally symmetric irreducible representations as well.

It is not a coincidence that the PLMs belong to the totally symmetric irreducible representation of the isolated donor and acceptor moieties. This connection can be understood by analyzing the corresponding changes in the electron density that accompany the transition. For example, the PLM on the naphthalene part shown in Fig. 2 largely corresponds to a symmetric ring-stretching mode, involving carbons C1-C2, C3-C4, C5-C6, and C7-C8. During the exciton transfer, an electron is promoted from the naphthalene HOMO to the naphthalene LUMO (shown in Fig. 5) and one expects that changes in the bond lengths should reflect the changes in electronic population. For naphthalene, a HOMO → LUMO transition would decrease the π-bond order between C1-C2, C3-C4, C5-C6, and C7-C8. Similar statements can be made for the phenyl and anthracene systems. In Sec. IV, we examine the PLMs and bond-length changes that occur during an energy transfer event. We show that, in a more quantitative way, while the PLMs can be understood in terms of geometric changes in the molecule, the reverse is not true. It is not at all straightforward to determine the PLMs by taking the difference between the initial and final geometries of the molecule.

We have established that the primary Lanzcos mode is more like the normal modes of the isolated donor and acceptor moieties, rather than of the entire molecule. The relation between the PLM and bond length change partially verifies our statement. In this section, we explore the localization in greater detail. First, we project the PLM onto the normal modes of entire molecule, which is shown in Fig. 6, as a comparison to Fig. 3, where the normal modes of the individual moieties are used. It is clear that on the basis of entire molecule, more modes are significantly involved. This is in contrast to the projection of the PLM onto only the local modes of the isolated moieties, only a handful modes contribute. Since the PLM best captures the initial nuclear dynamics, its close relation with local sites reveals the fact that the exciton transfer here is more like a local excitation/de-excitation coupled by exchange interaction, in accordance with the Dexter mechanism.¹⁶ 

One may argue that bridge is not included in our comparison. As showed in Fig. 2, one can see that very few motions involving the bridge contribute to the PLM. Thus, we anticipate that we can construct the PLM using only local modes on the donor and acceptor units. In Sec. V, we quantitatively verify this by using only dominant modes on the donor and acceptor to construct coupling modes.

However, before we can compare different ways to construct the Lanczos modes, we first need to address a subtle problem. In Sec. III, we symmetrized the geometry of naphthalene for the convenience in assigning irreducible representations to the vibrational modes. However, for constructing the modes for the entire unit one needs to use the optimized rather than the symmetrized geometry. In Fig. 7, we show the project of the PLM onto the modes for optimized naphthalene. The irreducible representations have been roughly assigned to the dominant modes. The projections are almost same to the ones on the symmetrized geometry, except a new B_2u mode that is involved in c14ee. However, it does not exist in d26ae. The reason is that the most significant modes of c14ee and d26ae are slightly different. The corresponding displacement vectors are embedded in Fig. 7. For c14ee, the dominant mode is not an exact A_g mode, because carbons move in the way of A_g, while hydrogens behave like B_2u. On the contrary, d26ae has a mode more similar to A_g. As a result, in the projection of naphthalene in c14ee, a B_2u mode is needed to correct the hydrogen motions.

To construct coupling modes, we collect a small number of modes on both naphthalene and benzaldehyde in the order of decreasing weights, then sum up the displacement vectors multiplied by their weights. We then benchmark the quality of taking various different combinations of normal modes by computing the coupling auto-correlation function in Eq. (3) and comparing to the exact result obtained by including all modes. We denote these as NxBy, where x denotes the number of naphthalene local modes and y denotes the number of benzaldehyde local modes used in each case. NfBf denotes that all local modes are included. The correlation functions and rate constants derived from various combinations of normal modes are summarized in Fig. 8.

First, if we take all local modes from both moieties, (NfBf), the resulting approximated correlation function is indistinguishable from the exact correlation function. This verifies our observation that the bridge unit is not needed in constructing the PLM. In order to achieve adequate agreement with the exact correlation function, at least 10 of the strongest contributing local modes are needed from both the donor and acceptor units.

However, as seen on our previous work,¹² the approximate correlation function does not need to give perfect agreement with the exact correlation function to produce an accurate transition rate constant (Eq. (4)). To assess the accuracy of using different combinations of local modes, we calculate the transition rate constants for various combinations up to 10 modes on each moiety and compare to the exact rate. These data are summarized in the form of contour plots in Figs. 8(c) and 8(d). For the c14ee case, there is an optimal “valley” whereby taking 3 modes from benzaldehyde and 4-6 modes from naphthalene produces an agreement with the exact result to within 5%. For the d26ae case, fewer number of benzaldehyde modes than naphthalene modes are needed to give agreement to within 5% of the exact rate constant. This is easily anticipated since naphthalene is simply a larger molecule with more normal modes.

We showed that the PLMs agree with bond length changes on the moieties. We next consider the relation between the PLM and geometry change accompanying the transition. First, we use the PLM to approximate geometry change by distorting the molecule along the PLM. The magnitude of the distortion is determined by minimizing the average bond length difference of the four bonds on naphthalene which stretch significantly in PLM between the optimized and approximate geometries. In Table I, we list the differences between two geometries. Comparing the maximum difference of bond lengths and angles between the approximated and true final states, the differences are very small, so PLMs seem to give fairly good approximations of geometry changes.

However, this is only one side of the story and we need to try to do this in reverse, i.e., can we determine the PLM from simply geometric changes accompanying the transition? To see if this is possible, we can project the geometry change onto the normal modes of molecule and moieties, respectively. The results for our two test cases are shown in Figs. 9 and 10. When projected onto the whole molecule, we see that geometry change is mostly dominated by one very low frequency mode corresponding to an internal torsion of the entire donor-bridge-acceptor molecule. Compared to the PLMs in the same basis in Fig. 6, we see that PLMs have many more modes involved. However, if we compare the projections onto the moieties in Figs. 3 and 9, the geometry changes resemble PLMs well, except in the low frequency region. In essence, simply changing representation completely changes the implied dynamics.

The resolution to this inconsistency is to distinguish between the two types of geometric changes. One is the internal changes within the moieties. These are well represented by the PLM. The other is the gross motion of the entire molecule, which we shall refer to as the “RM.” This latter motion is not well represented by the PLM. These two kinds of motion correspond to two different steps that accompany the electronic transition. The internal motion primarily affects C=C bond lengths and hence is strongly coupled to the π-electronic degrees of freedom, while gross motion of moieties corresponds to relaxation after exciton transfer has occurred. This is verified in Table I, where we use the relaxation mode to approximate the final geometry. The procedure is similar to that of PLM. Although RM is the dominant mode in the projection, it gives much worse approximated bond lengths and angles.

Since the PLM explains the short-time dynamics and can approximate the local geometry changes within the individual moieties very well, we anticipate that the PLMs facilitate the energy transfer step. To verify our prediction, we computed correlation functions of the electron/phonon coupling operators (Eq. (4)) using different subsets of nuclear motions to estimate the couplings: PLM, RM, GC, and GC-RM component eliminated. The latter two correspond to taking the total geometric difference between the initial and final states of the molecule and projecting this onto normal modes and to taking the geometric difference and subtracting off the relaxation mode. The various correlation functions and rate constants are summarized in Fig. 11. The exact results are obtained by using all normal modes and couplings. The PLM gives nearly perfect agreement to the exact calculation for both the time-correlation functions (Figs. 11(a) and 11(b)) and the computed rates for both molecules considered (Fig. 11(c)). The RM gives no contribution to the correlation function or the rate. Finally, using just the over-all GC gives a poor estimate of both the correlation function and the rate. GC-RM gives the initial transient dynamics in the correlation function and a reasonable approximation to the overall rate. However, the best agreement overall is obtained by using the PLM identified using the Lanczos search algorithm.

We presented here a further exploration of the classes of nuclear motions that accompany the electronic energy transition. We show that deeper insight can be gained by analyzing the primary modes and comparing them to both the geometric distortion of the molecule and to the electronic orbitals involved in the transition.

Our results also suggest that the energy transfer process in these systems occurs in two distinct steps. The initial transfer from the donor to acceptor occurs in a fixed nuclear frame–as per the Condon principle. This fast transfer is facilitated by the PLM and involves excitations in the C=C bond lengths of donor and acceptor fragments in accordance with changes in the HOMO/LUMO populations of each fragment. After the exciton transfer, the entire molecule undergoes a slow but large geometric relaxation under the new electronic distribution.

Finally, we emphasize that the PLM approach and this analysis need not be limited to the energy transfer cases studied here. Combined with the time-convolutionless master equation approach, this provides a robust and efficient way for computing state-to-state transitions rates for a wide class of molecular systems.

This work was funded in part by grants from the National Science Foundation (No. CHE-1362006) and the Robert A. Welch Foundation (E-1337). We thank Dr. Hao Li for useful comments and suggestions.