Assessment of simple and intuitive semiempirical approximations for non-adiabatic coupling vectors in the frame of (LC)-TDDFTB

TL;DR

This paper introduces two semiempirical methods for approximating non-adiabatic coupling vectors in quantum chemistry, benchmarking their qualitative performance and exploring their implications for electronic coupling and fluorescence in extended molecular systems.

Contribution

It presents two simple semiempirical schemes for calculating non-adiabatic coupling vectors and evaluates their qualitative effectiveness within the (LC)-TDDFTB framework.

Findings

Atomic NAC vectors are often qualitatively reproduced.

Quantitative accuracy of NAC magnitudes is limited.

Size-dependent trends in fluorescence yields are qualitatively predicted.

Abstract

In this study, two different simple and intuitive semiempirical schemes for computing approximate non-adiabatic coupling vectors (NACVs) between the ground and excited electronic states are presented. The first approximation makes use of Mulliken transition charges, while the second is based on derivative coupling vectors between localized molecular orbitals. Both approximations lend themselves easily to implementation within a whole spectrum of semiempirical quantum-chemical semiempirical methods. Here we present the implementation within the tight-binding DFT and benchmark its performance against analytical TD-DFT NAC vectors for a range of planar fluorescent chromophores at the Franck-Condon point. The pattern of the atomic NAC vectors is often reproduced, but the relative magnitude and total length of the NAC vector are often in serious error. Although quantitative predictions are not possible, these simple and intuitive approximations allow to explain, in a qualitative way, trends in the electronic coupling in extended molecular systems and complex materials. In this context, we investigate how the non-adiabatic coupling depends on the delocalization length of an excitation in chromophoric oligomers based on a simple model. Finally, we make general qualitative predictions on the size dependence of the fluorescence quantum yields in extended molecular systems, and illustrate those on the example of triply fused porphyrin tapes with increasing length.

Tables

Table 1. Table 1 : Dependence of the vertical excitation energy E vert subscript 𝐸 vert E_{\text{vert}} , radiative rate k rad subscript 𝑘 rad k_{\text{rad}} , non-radiative rate k nr subscript 𝑘 nr k_{\text{nr}} and fluorescence quantum yield Q ​ Y = k rad k rad + k nr 𝑄 𝑌 subscript 𝑘 rad subscript 𝑘 rad subscript 𝑘 nr QY=\frac{k_{\text{rad}}}{k_{\text{rad}}+k_{\text{nr}}} on the number of porphyrin units n 𝑛 n in the triply-fused porphyrin tapes T n .

Tn	$E_{\text{vert}}$ / eV	$k_{\text{rad}}$ (s-1)	$k_{\text{nr}}$ (s-1)	$QY$
1	2.35	9.6e+05	3.0e-07	1.0e+00
2	1.51	1.5e+07	8.2e+07	1.6e-01
3	1.19	3.4e+07	3.2e+09	1.0e-02

Equations

$${\mathbf{\tau }}_{mn}=⟨{\Psi }_m\mid \frac{\partial {\Psi }_n}{\partial \mathbf{R}}⟩.$$

$$⟨{\Psi }_m\mid \frac{\partial {\Psi }_n}{\partial \mathbf{R}}⟩=\frac{⟨{\Psi }_m\mid \frac{\partial \hat{H}}{\partial \mathbf{R}}\mid {\Psi }_n⟩}{E_n-E_m}$$

$$⟨{\Phi }_i\mid \hat{H}\mid {\Psi }_n⟩={\tilde{E}}_n⟨{\Phi }_i\mid {\Psi }_n⟩\forall i=1,\dots ,N_{\text{basis}}$$

$${\hat{V}}_{ne}=\sum _A^{\text{atoms}}\sum _i^{\text{electrons}}\frac{-Z_A}{|{\mathbf{R}}_A-{\mathbf{r}}_i|},$$

$${\mathbf{\tau }}_{mn}^A=\frac{Z_A}{E_n-E_m}\int d\mathbf{r}\frac{{\mathbf{R}}_A-\mathbf{r}}{|{\mathbf{R}}_A-\mathbf{r}{|}^3}{\rho }_{mn}(\mathbf{r}),$$

$${\rho }_{mn}(\mathbf{r})=N\int \dots \int d{\mathbf{r}}_2\dots d{\mathbf{r}}_N{\Psi }_m^{\ast }(\mathbf{r},{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N){\Psi }_n(\mathbf{r},{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N).$$

$${\mathbf{\tau }}_{mn}^A=\frac{-Z_A}{E_n-E_m}\int d\mathbf{r}\frac{\nabla {\rho }_{mn}(\mathbf{r})}{|{\mathbf{R}}_A-\mathbf{r}|}.$$

$${\rho }_{mn}(\mathbf{r})\approx \sum _B{q}_B\delta (\mathbf{r}-{\mathbf{R}}_B).$$

$${\mathbf{\tau }}_{mn}^A\approx \frac{Z_A}{E_n-E_m}\sum _{B\ne A}{q}_B\frac{{\mathbf{R}}_A-{\mathbf{R}}_B}{|{\mathbf{R}}_A-{\mathbf{R}}_B{|}^3}.$$

$${\mathbf{\mu }}_{mn}\approx \sum _A{q}_A{\mathbf{R}}_A.$$

$$\mid {\psi }_i⟩=\sum _{\mu }{c}_{\mu i}\mid \mu ⟩$$

$$\sum _{\nu }\left(H_{\mu \nu }^0-{ϵ}_i{S}_{\mu \nu }\right)c_{\nu i}=0$$

$$H_{\mu \nu }^0=⟨\mu \mid H^{\text{KS}}[{\rho }_0]\mid \nu ⟩$$

$$S_{\mu \nu }=⟨\mu \mid \nu ⟩$$

$$\begin{array}{rl}{\mathbf{d}}_{ij}^A& =⟨{\psi }_i\mid \frac{\partial {\psi }_j}{\partial {\mathbf{R}}_A}⟩\\ & \approx \frac{1}{{ϵ}_i-{ϵ}_j}\sum _{\mu ,\nu }{c}_{\mu i}^{\ast }{c}_{\nu j}\left[-\frac{\partial H_{\mu \nu }^0}{\partial {\mathbf{R}}_A}+\frac{{ϵ}_i+{ϵ}_j}{2}\frac{\partial S_{\mu \nu }}{\partial {\mathbf{R}}_A}\right]\end{array}$$

$$\mid {\Psi }_n⟩=\sum _{o\in \text{occ}}\sum _{v\in \text{virt}}{C}_{ov}^{(n)}\mid {\Psi }_o^v⟩$$

$${\mathbf{\tau }}_{0n}^A=⟨{\Psi }_0\mid \frac{\partial {\Psi }_n}{\partial {\mathbf{R}}_A}⟩=\sum _{o\in \text{occ}}\sum _{v\in \text{virt}}{C}_{ov}^{(n)}{\mathbf{d}}_{ov}^A$$

$$⟨{\Psi }_m\mid \frac{\partial {\Psi }_n}{\partial \mathbf{R}}⟩\cdot \frac{d\mathbf{R}}{dt}\approx \frac{1}{\Delta t}⟨{\Psi }_m(t)\mid {\Psi }_n(t+\Delta t)⟩$$

$$k_{\text{rad}}$$

$$k_{\text{IC}}$$

$$\text{QY}=\frac{1}{1+\frac{k_{\text{IC}}}{k_{\text{rad}}}}$$

$$x_i$$

$$q_i$$

$$|\mathbf{\mu }|=\sum _{i=0}^{2M-1}{q}_i{x}_i=qh$$

$${\tau }^{(i)}=\frac{q}{M{h}^2}\sum _{j\ne i}^{2M-1}(-1{)}^j\frac{i-j}{|i-j{|}^3}.$$

$$|\mathbf{\tau }|=\sqrt{\sum _{i=0}^{2M-1}({\tau }^{(i)}{)}^2}.$$

$$V_{ne}(\mathbf{r})=\sum _A\frac{-Z_A}{|{\mathbf{R}}_A-\mathbf{r}|}.$$

$$\begin{array}{rl}{\mathbf{\tau }}_{mn}^A& =\frac{1}{E_n-E_m}\int d\mathbf{r}\left({\nabla }_A{V}_{en}(\mathbf{r})\right){\rho }_{mn}(\mathbf{r})\\ & =\frac{-Z_A}{E_n-E_m}\int d\mathbf{r}\left({\nabla }_A\frac{1}{|{\mathbf{R}}_A-\mathbf{r}|}\right){\rho }_{mn}(\mathbf{r})\\ & =\frac{-Z_A}{E_n-E_m}\int d\mathbf{r}\left((-\nabla )\frac{1}{|{\mathbf{R}}_A-\mathbf{r}|}\right){\rho }_{mn}(\mathbf{r})\end{array}$$

$$\begin{array}{rl}{\mathbf{\tau }}_{mn}^A& =\frac{Z_A}{E_n-E_m}\left\{\cancel{{\left[\frac{{\rho }_{mn}(\mathbf{r})}{|{\mathbf{R}}_A-\mathbf{r}|}\right]}_{\infty }}-\int d\mathbf{r}\frac{\nabla {\rho }_{mn}(\mathbf{r})}{|{\mathbf{R}}_A-\mathbf{r}|}\right\}\\ & =-\frac{Z_A}{E_n-E_m}\int d\mathbf{r}\frac{\nabla {\rho }_{mn}(\mathbf{r})}{|{\mathbf{R}}_A-\mathbf{r}|}\end{array}$$

$$\begin{array}{rl}|\mathbf{\mu }|& =\sum _{i=0}^{2M-1}{q}_i{x}_i=\sum _i(-1{)}^i\frac{q}{M}ih=qh\frac{1}{M}\sum _{i=0}^{2M-1}(-1{)}^{i}i\\ & =qh\frac{1}{M}\left\{\sum _{j=0}^M(-1{)}^{2j}2j+\sum _{j=0}^{M-1}(-1{)}^{2j+1}(2j+1)\right\}\\ & =qh\frac{1}{M}\left\{\sum _{j=0}^{M-1}2j-\sum _{j=0}^{M-1}2j+2M-\sum _{j=0}^{M-1}1\right\}=qh\end{array}$$

Full text

I Introduction

Non-adiabatic coupling vectors (NACVs) play an important role in photochemistry. They describe the coupling between Born-Oppenheimer surfaces due to the nuclear kinetic energy and allow transitions between electronic states in the absence of radiation Baer (2006). They are a vital ingredient in

• •

non-adiabatic molecular dynamics simulations (surface hopping on “on the fly”) Tapavicza et al. (2013),

• •

searching for minimal energy conical intersections Ragazos et al. (1992)

• •

and predicting non-radiative transition rates and fluorescence quantum yields Valiev et al. (2018).

The brute force method for computing non-adiabatic coupling vectors is numerical differentiation of the wavefunction with respect to the atomic positions, which requires at least $3N_{\text{atoms}}$ electronic structure calculations. In the context of TD-DFT exact coupling vectors can be obtained analytically Send and Furche (2010) in an efficient way, but the implementation of the method is complicated. Having a simple and intuitive approximation for NACVs that may be combined with any semiempirical electronic structure method is therefore highly desirable.

In this article we compare two different semiempirical approximation for calculating non-adiabatic coupling vectors between the ground state and an excited state (usually $S_{1}$), which have been implemented in the frame of tight-binding DFT Humeniuk and Mitrić (2017). The first approximation is based on transition charges: In analogy with the transition dipole moment, the NACV is obtained simply from the transition charges and the molecular geometry. The second approximation, which has been propounded by Abad et.al. Abad et al. (2013), is based on molecular orbitals: Non-adiabatic couplings between Kohn-Sham orbitals are constructed from gradients of the overlap and Hamiltonian matrix between localized atomic orbitals, which are readily available in tight-binding DFT, since the same quantities are needed for evaluation of the energy gradient. Abad et.al. tested their approximation in the vicinity of conical intersections, where the magnitude of the NACVs is largely determined by the small energy gap. At these photochemical funnels the non-adiabatic coupling diverges and the transfer of population between electronic states is usually extremely fast. It remains to be investigate how well the approximation performs when the $S_{0}-S_{1}$ energy gap is large such as at the $S_{1}$ minimum or the Franck-Condon point where the length of the NACVs and their orientation relative to the normal modes determines the non-radiative transition rate.

The article is structured as follows: After a brief description of the approximations (in sections II.1.1 and II.1.2), we investigate how the non-adiabatic coupling depends on the delocalization length of an excitation in chromophoric oligomers (in section II.2). We then graphically compare the direction and magnitude of the approximate NACVs with their exact counterparts for a range of organic molecules with bright $\pi\pi^{*}$ excitations (in section III.1). Finally we make some qualitative predictions of fluorescence quantum yields in porphyrin tapes (section III.2).

II Theory

II.1 Semiempirical approximations

The first-order non-adiabatic coupling vector between two electronic Born-Oppenheimer states $m$ and $n$ is

[TABLE]

The coupling vector may be expressed as

[TABLE]

by differentiating the electronic Schrödinger equation on both sides with respect to the nuclear coordinates $\boldsymbol{R}$ and multiplying by $\langle\!\Psi_{m}\!\mid$ for $m\neq n$ and rearranging.

The derivation of this expression requires that $\hat{H}\mid\!\Psi_{n}\!\rangle=E_{n}\mid\!\Psi_{n}\!\rangle$ is satisfied exactly, which is a much stronger statement than just requiring that Schrödinger’s equation is satisfied after projecting onto a finite basis set ${\{\mid\!\Phi_{i}\!\rangle\}_{i=1,\ldots,N_{\text{basis}}}}$:

[TABLE]

Therefore eqn. (2) is strictly correct only if a complete basis set is used. In finite basis sets additional Pulay terms Pulay (1969) have to be considered which arise from the dependence of the basis set on the nuclear coordinates.

Nevertheless it is a good starting point for semiempirical approximations.

II.1.1 Approximation based on transition charges

Since the electronic Hamiltonian depends on the nuclear geometry only through the Coulomb attraction between nuclei and electrons,

[TABLE]

the coupling vector  (2) on atom $A$ simplifies to

[TABLE]

where we have also introduced the transition density matrix

[TABLE]

By partial integration of eqn. (5) (see appendix A) the NACV turns into

[TABLE]

This expression is very instructive since it shows that the coupling vector density is proportional to the gradient of the transition density. The largest contribution comes from points where $\boldsymbol{r}\approx\boldsymbol{R}_{A}$ due to the singularity of the Coulomb potential. Therefore we can say qualitatively that the non-adiabatic coupling vector on atom $A$ is approximately proportional to the gradient of the transition density around that atom.

To derive a semiempirical approximation for $\tau^{A}_{mn}$ let us return to eqn. (5) and assume that the transition density may be approximated by atomic transition charges (monopoles)

[TABLE]

where $\delta(\cdot)$ is Dirac’s $\delta$-function. This approximation is frequently employed in semiempirical methods such as tight-binding DFT Koskinen and Mäkinen (2009). The transition charges $q_{A}$ may be fitted to reproduce the electrostatic potential generated by the transition density (using the CHELPG algorithm) Madjet, Abdurahman, and Renger (2006) or they may be calculated as Mulliken transition charges from the transition density matrix. Substituting the monopole approximation  (8) into eqn.  (5) and using the property of the $\delta$-function, $\int\delta(x-x_{0})f(x)dx=f(x_{0})$, we get

[TABLE]

The term where $A=B$ was excluded to avoid dividing by zero. Only valence electrons are usually included in semiempirical calculations. Then the bare nuclear charge $Z_{A}$ should be replaced by the charge of the atomic core $Z_{A}^{\text{core}}$ (nucleus and core electrons), for instance in the case of carbon $Z_{A}^{\text{core}}=4$ instead of $Z_{A}=6$.

This approximation is completely analogous to how the transition dipole moment is calculated from the transition charges in the frame of TD-DFTB Niehaus et al. (2001),

[TABLE]

The simplicity of the derived approximate expressions enables us to make some general statements about the properties of the NACVs. The direction and length of NACVs can be deduced qualitatively by inspecting the transition density or the distribution of the transition charges:

• •

Coupling vectors are non-zero only on atoms which take part in an excitation.

• •

The coupling vectors point roughly along the direction where the transition density changes most strongly. Thus, if there is a node in the transition density between two atoms, the NACV on the atom is perpendicular to the nodal surface.

As a simple example consider the $\pi\pi^{*}$ excitation in ethene (Fig. 1). The transition charge is positive on one carbon, negative on the other and almost zero on the hydrogen atoms. Therefore the coupling vectors on the hydrogen atoms are zero. The transition charges change strongly from $+q$ to $-q$ when moving from one carbon to the other along the C$=$C bond. Therefore the NACVs on the carbons point along this bond.

The approximation fails completely when the transition density cannot be adequately described by monopoles. For instance in water, the HOMO-LUMO transition, $4a_{1}\leftarrow 1b_{1}$, has lobes of opposite sign below and above the molecular plane. The gradient of the transition density points perpendicularly to the molecular plane and is orthogonal to all vectors $\boldsymbol{R}_{A}-\boldsymbol{R}_{B}$. This implies that the coupling vector cannot be represented in the basis of bond vectors. The Mulliken transition charges are all zero, as is the approximate non-adiabatic coupling vector. In this case the approximation for the electric transition dipole given in Eq. 10 is also incorrect.

II.1.2 Approximation based on molecular orbitals

Here we briefly recapitulate how NACVs are calculated in the local-orbital scheme proposed in Ref. Abad et al. (2013) using the language of tight-binding DFT (DFTB). In DFTB a minimal basis set of valence atomic orbital is used. The molecular orbitals (MO) are linear combinations of these localized basis functions $\mid\!\mu\!\rangle$:

[TABLE]

The coefficients $c_{\mu i}$ for the molecular orbital $i$ are the eigenvector of the Kohn-Sham equation belonging to eigenenergy $\epsilon_{i}$:

[TABLE]

Matrix elements of the Kohn-Sham Hamiltonian at the reference density,

[TABLE]

the overlap matrix elements

[TABLE]

and their gradients are obtained from Slater-Koster rulesSlater and Koster (1954).

With the help of eqn.  (12) the authors of Ref. Abad et al. (2013) derived an approximate expression for non-adiabatic coupling vectors between molecular orbitals:

[TABLE]

In time-dependent density functional theory and its tight-binding version, excited states are represented as linear combinations of singly excited Slater determinants:

[TABLE]

Non-adiabatic coupling vectors between the many-electron ground and excited states are obtained by contraction of the single-particle coupling vectors $\boldsymbol{d}_{ij}^{A}$ with the coefficients $C^{(n)}$:

[TABLE]

It is worthwhile to highlight some of approximations made in the above derivation: (a) Eqn.  (12) neglects the dependence of the Hamiltonian on the density ($H[\rho]\approx H[\rho_{0}]$). This allowed to derive the relatively simple expression  (15) for the coupling vectors. However, in our calculations we use expression  (15) with MO coefficients obtained from solving the Kohn-Sham equations self-consistently. (b) In principle, the NACV contains a contribution from changes of the coefficients $\frac{\partial C^{(n)}_{ov}}{\partial\boldsymbol{R}}$, which is neglected. (c) The exact NACV diverges when the ground and excited state cross ($E_{1}\approx E_{0}$), whereas for the approximate NACV this happens when the HOMO-LUMO gap closes ($\epsilon_{\text{HOMO}}\approx\epsilon_{\text{LUMO}}$). It is well-known that HOMO-LUMO gaps are often significantly larger than $S_{0}-S_{1}$ excitation energies obtained from TD-DFT due to the mixing of single excitations. This suggests that the approximation will work best when such many-body effects are small so that the $S_{0}-S_{1}$ transition predominantly has HOMO $\to$ LUMO character.

Expression  (15) is ideally suited for tight-binding DFT, since the gradients of the matrix elements can be constructed very efficiently at runtime from precalculated Slater-Koster tables. Since the same quantities are needed for assembling the gradient of the energy, which is needed in any molecular dynamics (MD) simulation, the computation of the NACVs comes at little additional cost. This should be contrasted with the computational cost of an alternative method for computing the non-adiabatic couplings in MD simulations: In the surface hopping method Tully (1990) the electronic populations depend only on the scalar product between the NACV and the nuclear velocity vector. This scalar can be obtained directly from the overlap of the electronic wavefunctions at neighbouring timesteps Mitric et al. (2009) obviating the need for computing the NACVs:

[TABLE]

However, since each excited state is a linear combination of Slater determinants, the evaluation of the overlap entails a large number of determinants, rendering this scheme very expensive for large molecules, unless cutoff thresholds are used for culling determinants which contribute little to the overlap integral.

II.2 Qualitative fluorescence quantum yield

Approximations  (9) and  (10) provide qualitative guidelines on how to tune the electronic wavefunctions for increasing the fluorescence quantum yields. At the moment we will only focus on electronic effects, although vibrational effects can also be very important, as will become clear later on.

If the vibrational wavefunction is neglected, according to Fermi’s Golden rule the rates for radiative (spontaneous emission) and non-radiative (internal conversion) decays are proportional to the lengths squared of the transition dipole and non-adiabatic coupling vectors, respectively, between the ground state $S_{0}$ and the first excited state $S_{1}$:

[TABLE]

To increase the fluorescence quantum yield

[TABLE]

$k_{\text{rad}}$ needs to be maximized while $k_{\text{IC}}$ needs to be minimized.

This can be achieved by

• •

increasing the length of the transition dipole

and/or

• •

avoiding conical intersections, where $E_{1}=E_{0}$

• •

and reducing the gradient of the transition density.

To avoid the crossing of the energy levels of $S_{1}$ and $S_{0}$, the geometry should be rigid, so that we can assume there is a stable minimum on $S_{1}$ and the reorganization energy is small. Then it remains to maximize the length of the transition dipole moment and to minimize the gradient of the transition density. Since it is easier to analyze only one factor, we build a simple model, where the transition dipole is constant and only the length of the NAC vector changes.

II.2.1 1D model fluorophore

Consider a linear molecule (e. g. a polyene) with $2M$ atoms on an equidistant grid with spacing $h$ (see Fig. 2a). For simplicity, we assume that each atom has a single $p_{z}$ orbital and contributes one electron. The $S_{1}$ state is a HOMO-LUMO transition. The $S_{0}-S_{1}$ transition density has nodes between all atoms, so that the transition charges alternate between positive and negative values.

The atomic positions and transition charges for atom $i$ are given by

[TABLE]

For simplicity we set the nuclear charge to $Z=1$ and the excitation energy to $E_{1}-E_{0}=1$.

The transition dipole moment is independent of the number of atoms (see appendix B):

[TABLE]

The non-adiabatic coupling vector on the $i$-th atom is given by

[TABLE]

The length of the total non-adiabatic coupling vector is given by

[TABLE]

and needs to be evaluated numerically. If we plot the length of $\boldsymbol{\tau}$ against the number of atoms $2M$ that participate in the excitation (Fig. 2b), we see that the rate for internal conversion can be minimized by spreading the transition charge over as many atoms as possible while maintaining the same electric transition dipole moment. Since the transition charges change sign every second atom, the gradient of the transition density can be reduced only if the charges themselves are small. In order to keep the same transition dipole moment, the number of atoms over which the excitation is delocalized needs to be increased.

III Results

III.1 Comparison between approximate and exact NACVs

The two approximations for NACVs are tested for a series organic molecules with bright $\pi\pi^{*}$ transitions. Many of the selected molecules are fluorescent dyes which have a stable lowest excited singlet state (with the exception of the polyenes).

After optimizing the geometries at the AM1 level of theory, the lowest bright excited state was computed with TD-$\omega$B97XD/def2-SVP using Gaussian 16 Frisch et al. (2016). Analytical NACVs were obtained in the frame of TD-DFT Send and Furche (2010) via the keyword TD=NAC. These vectors serve as “exact” reference values against which the quality of the approximate vectors is measured. Approximate NACVs based on either Mulliken transition charges (according to eqn.  (9)) or localized orbitals (according to eqns.  (15) and  (17)) were computed in the frame of long-range corrected tight-binding DFT Humeniuk and Mitrić (2017). The comparison between the three types of NACVs is presented in a graphical way in Figs. 3 to 8 below. The components of the NAC vectors on each atom are shown as little red arrows. Since eigenvalue solvers produce eigenvectors with arbitrary global signs, only the relative orientation of the vectors to each other is important. A sign change in either the bra or the ket wavefunction is equivalent to flipping all vectors simultaneously. The NACVs were scaled by a factor (which is indicated in the upper left corner) so that the largest vectors in each figure has approximately the same length. We proceed by analyzing the quality of the semiempirical approximations as compared to the exact NACVs for each class of molecules:

Trans-polyenes (Fig. 3) are the simplest conjugated systems which behave like the linear 1D model discussed above. The lowest bright state ($B_{u}$) is polarized along the molecular axis. The transition density has nodes between neighbouring carbon atoms, which coincide with the positions where the tip of one arrow touches the tail of the next. The individual arrows become shorter as the number of carbon atoms increases from ethene to hexatriene, reflecting the smearing out of the transition charges over a larger area. The transition charge (TC) approximation overestimates the NACVs by a factor of 5 but gets the orientation of the vectors right. In turn, the localized orbital (LO) approximation underestimates the NACVs by at least a factor of 10 and fails to predict the orientation.

The cyanines Cy$N$ (Fig. 4) are fluorescent cationic dyes that consist of a polymethine chain connecting two nitrogens which are part of an indole moiety. Cy3, Cy5 and Cy7 differ by the number of carbon atoms in the bridge, in Cy3B Cooper et al. (2004) the polymethine chain is stabilized against deformation by additional aliphatic rings. In all cyanines the lowest bright excitation is localized on the polymethine chain, and consequently the NACVs are also limited to this region of the molecule. In the polymethine chain the orientation of the arrows alternates as is expected based on the location of the nodes in the transition density. The LO approximation predicts the position and orientation of the NAC vectors correctly but underestimates their magnitude by a factor of 3. The TC approximation yields NACVs that are spread out too much over non-chromophoric parts of the molecule, such as methyl groups in Cy3-Cy7 or the aliphatic rings in Cy3B. In the polymethine chain the NAC vectors all point in the same direction, but the total magnitude of the NACV is approximately correct.

Dicyanovinyl-substituted squaraines (Fig. 5) Mayerhöffer et al. (2013) are another class of fluorescent dyes. In squaraine-O position 3 of the indole moiety is replaced by oxygen, whereas in squaraine-CMe a methyl group is added. The excitation is localized on the central four-membered ring and the adjacent methine groups. The LO approximation reproduces the orientation of the vectors accurately, except for those on the C$\!=\!$O group, which are far too short. The TC approximation places the largest NAC vectors on two opposite carbon atoms in the four-membered ring, although the coupling vectors at these positions should be zero. As with the cyanines a tendency of TC is observed to place large NAC vectors on atoms that are not part of the chromophore.

Finally a selection of polycyclic aromatic hydrocarbons is considered in Figs. 6 and 7. The couplings were calculated for the lowest excited state. Since the ordering of states can be method-dependent, the symmetry label is given in brackets.

The ring systems give rise to complex patterns in the distribution of NAC vectors. In fluorene ($B_{2}$) and phenanthrene ($B_{2}$) the arrows are arranged in cycles around the outer six-membered rings. This pattern is reproduced by the LO approximation, whereas the TC pattern has no similarity with the exact results. In pyrene ($B_{1u}$), perylene ($B_{1u}$) and rubrene ($B_{1}$) and relative orientation of the vectors is reproduced correctly both by the LO and the TC approximations, however the relative magnitudes of the vectors differ considerably. The total magnitude of the coupling is severely underestimated by the LO approximation (by a factor of 3-10) and overestimated by the TC approximation (by as much as a factor of 10). In rubrene the excitation is strictly confined to the tetracene core. Inspite of this, the semiempirical approximations yield large vectors on the adjacent phenyl groups, which are perpendicular to the central tetracene.

III.2 Porphyrin tapes

We will now test the predictions of the 1D model from section II.2.1 for the porphyrin tapes that were synthesized by the Tsuda group Tsuda and Osuka (2001). These tapes consist of triply-fused zinc-porphyrins (the structure is shown as an inset in Fig. 9). The monomer units are linked through conjugation allowing the electrons to delocalize freely over the entire tape like particles in a box. The delocalization is reflected in the lowering of the excitation energy far into the infrared with increasing length. At the same time, delocalization of the transition density should also impact the magnitude of the electronic non-adiabatic coupling.

The transition dipole moments and NACVs were computed with long-range corrected TD-DFTB for the lowest $B_{1u}$ state, which is polarized along the long axis of the tape. For the monomer, dimer and trimer the NACVs are depicted in Fig. 8. As the conjugation extends over all porphyrin units, the transition dipole moment $\boldsymbol{\mu}$ grows approximately linearly with the size of the tape. However, since the transition charges are spread out over a larger area, the non-adiabatic coupling $\boldsymbol{\tau}$ grows sublinearly and saturates. The ratio between the lengths of the two vectors is shown in Fig. 9. Since the tapes are also very rigid, one can expect that ultrafast internal conversion through conical intersections, which usually requires some local deformation of the geometry, is not the dominant non-radiative decay channel. Based on this analysis one would expect the long tapes to have extremely high fluorescence quantum yields.

However, this is not the case. Article Cho et al. (2002) actually shows that the non-radiative rate increases rapidly with the length of the tapes so that the fluorescence is quenched as compared to the monomer. At first this appears to contradict the fact that the ratio of the electronic non-adiabatic coupling to the transition dipole moment, $\tau/\mu$, decreases. However, the reason for the fluorescence quenching in the oligomers is the lower energy gapTittelbach-Helmrich and Steer (1995) and the vastly higher density of states. The reduction in the electronic non-adiabatic coupling per porphyrin unit is more than compensated by the increase of accessible final vibrational states.

To verify this explanation we computed the radiative and non-radiative rates for the smallest porphyrin tapes using Fermi’s Golden Rule in the harmonic approximation following the steps of Ref. Valiev et al. (2018). 111Calculation of rates is based on the following approximations: The $S_{0}$ and $S_{1}$ potential energy surfaces are harmonic, share the same normal modes and frequencies but have different equilibrium geometries (shifted harmonic oscillators). Frequencies and normal modes are determined from a frequency calculation using $\omega$B97XD/dev2-SVG at the $S_{0}$ minimum. The Huang-Rhys factors are deduced from the gradient on $S_{1}$ at the Franck-Condon point. Total rates are obtained by summing over all transitions that start in the initial vibrational ground state on $S_{1}$ and lead to a final vibrational state on $S_{0}$. Radiative transitions may lead to any vibrational state lower in energy, while in non-radiative transitions the final vibrational states on $S_{0}$ have to be approximately isoenergetic with the initial vibrational state on $S_{1}$. Following Ref. Santoro et al. (2007) final states are grouped into classes Cn depending on the number $n$ of simultaneously excited modes. The summation is limited to classes $C_{0}-C_{8}$, which captures most of the radiative rate and a large fraction of the non-radiative rate. Modes are sorted in decreasing order by Franck-Condon factors and the number of modes from which excitations are allowed, is reduced, until there are no more than $10^{8}$ elements per class left. For the smallest tapes optimizations and frequency calculations on the ground state and the first excited state with $1B_{1u}$ symmetry at the TD-DFT level of theory are still feasible. The resulting rates and quantum yields are listed in table 1:

The non-radiative rate jumps by orders of magnitude from the monomer to the dimer and increases further in the trimer. With the non-radiative rate increasing much faster than the radiative rate the quantum yield drops to zero, as observed in experiment.

Since the sum over final vibrational states necessarily has to be truncated, the reported non-radiative rates are only a lower limit. Even then it is clear that the non-adiabatic coupling between vibrational states and the sheer density of states is responsible for the fluorescence quenching.

IV Discussion

Judging the quality of the NAC vectors by visual inspection can be misleading since it suggests there is more agreement than there actually is. The symmetry of NAC vectors is related to the symmetry of the excited state. The relative orientation of the vectors in molecules with high symmetry, is therefore largely determined by the irreducible representation. In trans-butadiene ($C_{2h}$), for instance, only the relative orientation of two out of four vectors is not already fixed by symmetry. According to TD-DFT these two vectors not related by symmetry should be parallel, but the localized orbital method yields an antiparallel orientation (see Fig. 3). The orientation is thus entirely wrong and the magnitude is also wrong by a factor of 10.

The localized orbital method tends to underestimate the magnitude of the vectors: In ethene the vectors are too short by a factor of 60, in the cyanine dyes by a factor of 3 and in the porphyrin tapes by a factor of 8. The large error for a system as simple as ethene is surprising. Eqn.  (7) and Fig. 1 showed that the non-adiabatic coupling in the $\pi\pi^{*}$ state is due to the gradient of the transition density which points along the C-C bond. The transition charge approximation fares a little bit better in predicting the magnitude of the coupling, but it fails in predicting the distribution of the vectors: In the cyanines the excitation is strictly localized on the polyene bridge, but large vectors can be found on two adjacent methyl groups. An extreme example of this are the porphyrin tapes, where the largest vector is placed on the zinc atom, which does not take part in the excitation at all.

Comparison between the two approximations is hindered by the fact that one is derived from eqn.  (7) (gradient of transition density) and the other from eqn.  (1) (gradient of excited state wavefunction), but the two expressions are only equivalent in the basis set limit, and the minimal valence basis of DFTB is a long way off a complete basis set. The LO approximation considers terms that arise because basis functions are attached to the nuclei (Pulay terms) but neglects changes of the excitation coefficients ($\frac{\partial C^{(n)}_{ov}}{\partial\mathbf{R}}$). The TC approximation is independent of the basis set and thus cannot account for Pulay terms.

Some of the errors relative to TD-DFT might also be due to the tight-binding approximations: Semiempirical transition charges, excitation energies and molecular orbitals, which enter the expressions for the TC and LO approximations, differ from their ab initio counterparts. However, those sources of error are of minor importance. In fact, if we feed our TC approximation with transition charges that were fitted to reproduce the electrostatic potential of the TD-DFT transition density (using the PSPFFT libraryBudiardja and Cardall (2011) for solving the Poisson equation and the CHELPG algorithm Breneman and Wiberg (1990)), the resulting vectors are very similar to the tight-binding results. The valence basis set employed in DFTB is also not to blame. With a minimal STO-3G basis set the resulting TD-DFT NAC vectors are indistinguishable from the def2-SVP results.

V Conclusion

Two simple semiempirical approximations for non-adiabatic coupling vectors between excited singlet states and the ground state were implemented in the frame of (LC)-TDDFTB and compared with TD-DFT coupling vectors as benchmarks for a set of planar chromophores with bright $S_{1}$ states. The TC approximation is based on excitation energies, atom-centered transition charges and geometric information. In the LO approximation the coupling between many-body states is calculated from the coupling vectors between molecular orbitals.

While easy to implement and highly efficient, both approximations are not accurate enough to predict the absolute magnitude of the non-adiabatic coupling vector. In particular the LO approximation underestimates couplings by one order of magnitude. Nevertheless, the region in the molecule where the coupling is large can often be identified. For a series of fused porphyrin tapes the reduction in the electronic coupling per porphyrin unit can be explained by the increasing delocalization of the excitation. As a general rule, spreading transition charges over a larger area reduces the electronic non-adiabatic coupling. This however, does not imply that the fluorescence quantum yield may be increased simply by enlarging the delocalization length, since larger $\pi$-system also have larger nuclear non-adiabatic couplings due to the increased density of states.

The upshot is that quantitative NAC vectors cannot be obtained with these simple approximations. The implementation of analytical coupling vectors in the spirit of Ref. Send and Furche (2010) can in principle be adapted to tight-binding DFT in analogy to the analytic gradients Heringer et al. (2007) but will require a major effort. The LO approximation is a first step in that direction. Without going to these lengths, the TC approximation might be improved upon by including higher multipoles to represent the transition density more faithfully away from the molecular plane.

Acknowledgements

A.H. and R.M. acknowledge funding by the European Research Council (ERC) Consolidator Grant DYNAMO (Grant No. 646737).

Appendix A Detailed derivation of eqn.  (7)

The single-electron part of the electron-nuclear attraction is

[TABLE]

[TABLE]

By the partial integration rule $\int_{V}(\nabla f)g=[fg]_{\partial V}-\int_{V}f(\nabla g)$ this becomes

[TABLE]

Appendix B Transition dipole moment in the linear molecule

The following little calculation shows that the transition dipole moment is a constant, independently of the number of atoms $2M$:

[TABLE]