A rigorous nonorthogonal configuration interaction approach for the calculation of electronic couplings between diabatic states applied to singlet fission

For the design of efficient singlet fission chromophores, knowledge of the factors that govern the singlet fission rate is important. This rate is approximately proportional to the electronic coupling between the lowest (diabatic) spin singlet state that is populated following photoexcitation state and a so-called TT state. The latter state is characterised by two triplets, each localised on one of two neighbouring molecules, which are coupled into a singlet. Here, we show the applicability of a nonorthogonal configuration interaction approach for the calculation of this electronic coupling. The advantages of this rigorous approach are that (1) the coupling can be calculated directly, (2) it includes important correlation and orbital relaxation effects, and (3) it has a clear chemical interpretation in terms of molecular states. This approach is applied to calculate the electronic coupling for a biradicaloid molecule, viz. the bis(inner salt) of 2,5-dihydroxy-1,4-dimethyl-pyrazinium. The biradicaloid molecule is, based on the energetic criteria, a promising candidate for singlet fission. We show that the electronic coupling between the molecules is also sufficiently large for singlet fission, rendering molecules based on this chemical moiety interesting singlet fission chromophores. 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
In the development of solar cell technology, finding new materials that can convert solar energy into electricity with high efficiency is a very active area of research [1][2][3][4]. Materials that can generate multiple electron-hole pairs per single absorbed photon have potential to improve the efficiency of solar cells [4,5]. This process is known as multiple exciton generation. An example of multiple exciton generation is singlet fission (SF), which enables the exploitation of high energy photons with a minimum of thermal energy loss [6,7]. It is a fast (sub-ps) and radiationless process [8], in which the first excited singlet state of a chromophore (S 1 ) transfers part of its energy to the neighbouring ground state chromophore (S 0 ) resulting in two (local) triplets that are coupled into a singlet ( 1 TT state, Fig. 1). This is a spin allowed process, as the total spin S = 0 of the S 0 S 1 and 1 TT states is conserved upon transition from the S 0 S 1 state to the 1 TT state [8].
For an efficient SF process a number of criteria have to be fulfilled. First, the SF process has to be faster than other competing processes. Ideally this process should be isoergic or slightly exoergic E(S 1 ) ! 2E(T 1 ). There are some guidelines that have been studied and applied for designing SF chromophores that fulfil this criterion [9]. Second, for a fast SF process, there should be some interaction between the two chromophores, but it should not lead to charge transfer state or excimer formation [9]. Meeting these criteria in practice is difficult, but it is not impossible to investigate them and study the SF mechanism rigorously with the use of theoretical chemistry and computational modelling.
Giving a full (and correct) computational description of the SF process is challenging. For example, in many molecules that could show SF, the excited S 1 state has double excitation character. Density Functional Theory (DFT) and Time-Dependent DFT can be used to determine the equilibrium geometry of this state, but these methods are not suitable for determining its relative energy. However, multiconfigurational methods such as Complete Active Space Self Consistent Field (CASSCF) and Restricted Active Space SCF (RASSCF) followed by second order perturbation theory are suitable to determine the energies, but are computationally too demanding for the geometry optimization. h jhW i jHjW f ij 2 qðEÞ [8], in which the wavefunctions of the initial and final diabatic states are simply those of the S 0 S 1 and 1 TT states. In this approximation, the electron transfer rate is determined by the electronic coupling between the diabatic states and the density of states factor. The present study will focus only on the calculation of the electronic coupling matrix element hW i jHjW f i. This matrix element can be evaluated in various ways, using cf. a ZINDO/CISD approximation [10]. Alternatively, one can use phenomenological models in which only the frontier molecular orbitals of interacting chromophores are considered [8], or DFT [11]. An approach based on localisation of the frontier molecular orbitals followed by transformation of the Fock matrix to this basis to determine the coupling has also been used to study the effects of vibrations on the electronic coupling in covalent tetracene dimers [12]. Recently, a nonorthogonal model to calculate the geometry dependence of the coupling between ethene dimers has been introduced, which considers the HOMO and LUMO only [13]. However, the existing approaches lack accuracy. These models ignore orbital relaxation and/or do not include electron correlation in a systematic way.
In this study we introduce a rigorous nonorthogonal configuration interaction approach that enables us to calculate the electronic coupling matrix element directly, using state specific CASSCF wavefunctions for each of the chromophores involved in the SF process. Our approach is able to incorporate important electron correlation and orbital relaxation effects. It is based on the use of a scheme introduced earlier [14,15], in which we express the diabatic states of an ensemble of molecules in terms of nonorthogonal configuration interaction (NOCI) wavefunctions. The manyelectron basis states in which these NOCI wavefunctions are expanded are formed by antisymmetrized products of molecular wavefunctions of the CASSCF type. In a previous study [14], due to technical limitations, such a basis state was approximated by computing the CASCI wavefunction of an ensemble of molecules with the orbitals obtained from CASSCF calculations on the individual molecules. The CASCI wavefunction contains unwanted, contaminating, charge transfer contributions that thwart the interpretation of the basis state as being composed of several molecular states. Moreover, these charge transfer contributions lead to long CI expansions, especially when many molecules are included in an ensemble. Furthermore, the CASCI calculation requires an unwanted intermediate (Löwdin) orthogonalization of molecular orbitals, which also obscures the interpretation of the state. Contrary to this previous approximate ansatz, in the present study the many-electron basis states in the NOCI wavefunc-tions are proper antisymmetrized products of molecular wavefunctions of the CASSCF type. The orbitals in these basis states are the molecular orbitals, without any orthogonalizations, as the CASCI step to obtain the CI coefficients is no longer needed.
The main advantages of this rigorous implementation compared to previous models, using orthogonal approaches, for calculating the electronic coupling between diabatic states are: (i) the explicit computation of the Hamiltonian matrix elements, (ii) the systematic inclusion of non-dynamical electron correlation and orbital relaxation effects, (iii) a clear chemical interpretation of the states involved, and (iv) compactness of the wavefunction.
We take an initial state i j i of an ensemble of neighbouring molecules as IJKL j i where I, J, K and L indicate the ground state wavefunction of the molecules in the ensemble. This ensemble will in the following be denoted as ''cluster". An intermediate, local where J T and K T are excited triplet states on neighbouring molecules J and K, which are coupled into a singlet.
The study of the delocalization of the excited singlet state involves computation of hm J jHjm K i and the study of the transition rate between m J and f j i involves the computation of hm J jHjf i where H is the Hamiltonian of the cluster. The computation of these matrix elements is non-trivial because the orbitals of the molecular wavefunctions in m J , m K j i, and f j i are different and mutually nonorthogonal. An interesting question that can be addressed with this approach is whether charge transfer states such as IJ þ K À L play a role in the SF process, either as intermediate or as virtual states. Electronic relaxation effects are very important in the excitation, delocalization and fission processes and the orbitals that are optimal for I, for I S , for I T , for I + and for I À are all quite different. It is therefore preferable both from a computational and from a conceptual viewpoint, to express each molecular state, and therewith each of these diabatic states, in its own optimised orbital set.
In our approach, basis states of the types and IJ þ K À L are used as the many electron basis for the diabatic states. These cluster basis states describe the ground states of all molecules in the cluster, singlet excitations on one of the molecules, triplet excitations on two neighbouring molecules, and intermolecular charge transfer between molecules. The cluster basis states are constructed as antisymmetrized products of the molecular wavefunctions, with optimised orbitals for each state of each molecule. This implies that the orbitals of different molecular wavefunctions in the cluster basis states are different and mutually not orthogonal. Hence, the construction of the cluster basis states and computation of the matrix elements are non-trivial.
To illustrate the applicability of this approach, we discuss the calculation of the electronic coupling matrix element between the nonorthogonal diabatic S 0 S 1 and 1 TT states for a molecule proposed by Michl and co-workers as a potential SF chromophore [16]. This so-called biradicaloid molecule (Fig. 2) has been selected and proposed to be synthesized because, based on quantum chemical calculations, it has been found to fulfil the excitation energy criterion E(S 1 ) % 2E(T 1 ) for a potential SF chromophore [16]. In addition, we show for this system the effect of the arrangement of neighbouring chromophores on the computed electronic coupling. In our previous work on tetracene [14], it was shown that only the nearest neighbour couplings are significant, and that these are insensitive to the cluster size. In addition, an efficient SF was observed in solution of one photoexcited and one ground state of TIPS-pentacene, showing that the involvement of two chromophores is sufficient to detect the SF [17]. Therefore, for the present purpose it suffices to use small clusters consisting of only two neighbouring molecules.

Methods
The crystal structure of the biradicaloid molecule is unknown, therefore, to determine a possible crystal structure, periodic DFT calculations starting from the known crystal structure of a related compound, namely 2,5-dimethyl-1,4-benzoquinone [18], were performed using the CRYSTAL14 code [19]. The PBE functional and the 6-21G basis set were employed. In addition, the Grimme dispersion correction [20] was included for the optimization of atom positions and cell parameters.
The excitation energies for the biradicaloid molecule were calculated using two different active spaces, i.e. CASSCF(2,2) and CASSCF (6,5). Dynamical correlation was included up to the second order perturbation theory (CASPT2) following the CASSCF(6,5) calculation. All CASSCF and CASPT2 calculations were performed using the ANO-L basis set contracted to C,N,O[3s2p1d]/H[2s] [21]. All of these calculations were performed using the MOLCAS 7.4 code [22].
For technical reasons, the CASSCF(2,2) wavefunctions of the ground state, first excited singlet and triplet states, as well as the ROHF wavefunctions of the cation and anion were also computed using the GAMESS-UK code [23]. These wavefunctions were subsequently used to construct the cluster basis states. The clusters consist of two molecules whose wavefunctions A and B can be combined to form six different antisymmetrized product wavefunctions, which are the cluster basis states. They are listed as follows: one basis state describing the ground states on both molecules W S 0 S 0 ¼Â½AB, two basis states describing the singlet excitation on one of the molecules W S 0 S 1 ¼Â½AB S and W S 1 S 0 ¼Â½A S B, one basis state describing two triplets that are coupled into a singlet W1 TT ¼Â½A T B T , and two basis states describing These cluster basis states were used to determine the diabatic S[1] and S [2] states and the 1 TT state. The Hamiltonian and overlap matrix elements between these nonorthogonal cluster basis states were calculated using the GNOME code that had been developed earlier in our group [24].
First, diabatic S [1] and S [2] states were obtained from a 2 Â 2 NOCI calculation in the basis of two cluster basis states, each having one molecule in its first excited singlet state. The diabatic 1 TT state is the 1 TT basis state. To investigate the effect of charge transfer states on the diabatic states, the charge-transfer basis states were added to the NOCI calculation for the S[1] and S [2] states (giving a 4 Â 4 NOCI) and to the NOCI calculation for the diabatic 1 TT state (giving a 3 Â 3 NOCI). To investigate the importance of each cluster basis state in the diabatic S [1], S[2] and 1 TT states, the weights (W i ) of the cluster basis states i, were calculated using the Gallup and Norbeck scheme [25], W i ¼ jc i j 2 =ðS À1 Þ ii , where c i is the CI coefficient of the basis state i and ðS À1 Þ ii is the ii th element of the inverse of the overlap matrix.

Results and discussion
The periodic DFT calculations of the biradicaloid showed a P À1 symmetry with the absence of imaginary frequencies and provided the final cell parameters: a = 3.578 Å, b = 8.757 Å, c = 9.413 Å, a = 96.55°, b = 96.77°, c = 109.11°, and q = 1.703 g/cm 3 . Fig. 3 shows the resulting crystal structure. Since to the best of our knowledge the molecule has not been synthesized yet, there are no experimental data to compare with. There are two different stacks in the crystal structure, to be denoted stack A and stack B. The band structure (not shown) shows dispersion mainly in the stack directions, thus we considered two intra-stack pairs of neighbouring molecules (in stack A and stack B, respectively). We also considered an A-B inter-stack pair. The intra-stack pairs show plike stacking (slip-stack) of two molecules while the inter-stack pair shows an arrangement of two neighbouring molecules, one is taken from stack A and the other one is taken from stack B. One difference between stack A and stack B is the N-N distance between two molecules, i.e. 3.854 Å and 3.638 Å, respectively (see Fig. 3).
An interesting property of this biradicaloid molecule is the multireference character of the ground state. A CASSCF(2,2) calculation gives natural orbital occupation numbers of 1.76 and 0.24 for the p-type HOMO and the p ⁄ -type LUMO, respectively, and hence, single reference methods are not suitable even to describe its ground state. The first excited singlet and triplet states of the biradicaloid molecule have mainly a single excitation from HOMO to LUMO and  therefore have p ? p ⁄ character as can be seen in Fig. 4. The excitation energies of the first excited singlet and triplet states computed with CASSCF(2,2) and CASSCF(6,5)/CASPT2 are listed in Table 1. The CASSCF(6,5)/CASPT2 excitation energies for the first excited singlet state is higher than the value reported by Akdag, et al. [16], in which a larger active space had been used [CASSCF (22,14)], in contrast with the CASSCF excitation energy for the first excited triplet state. The inclusion of dynamical correlation by perturbation theory up to the second order (CASPT2) improves the excitation energies considerably. The excitation energy of the first excited singlet state computed with CASSCF(2,2) is 1.0 eV higher than those computed with a larger active space. Since an important objective of this study is to prove the principle of the approach, and we are interested in estimating the order of magnitude of the electronic coupling, we used the CASSCF(2,2) wavefunctions to construct the cluster basis states for the diabatic S [1] and S [2] states and the 1 TT states even if their relative energies are rather poor.
As described previously, the combination of different molecular wavefunctions by taking their antisymmetrized product, gives six cluster basis states that were used to determine the diabatic S [1] and S [2] states and the 1 TT state. The relative energies of these cluster basis states for the pairs taken from stack A and stack B, and the inter-stack pair are listed in Table 2. The cluster basis states W S 0 S 1 and W S 1 S 0 describe the singlet excitation on one of the molecules and were used as the basis to construct the diabatic S[1] and S [2] states. The relative energies of these cluster basis states are close to the excitation energy of the first singlet excited state of the molecule. The cluster basis state W1 TT represents by itself the diabatic 1 TT state. The relative energy of this cluster basis state is about twice the excitation energy of the first excited triplet state of the molecule. The two cluster basis states W CT1 and W CT2 describe the charge transfer states in which an electron is transferred from one molecule to another. In all cases the relative energies of these charge-transfer basis states are 4.5 eV or more above the ground state energy. In principle these charge-transfer basis states may play a role in the SF process either as intermediate or as virtual states, but due to their high energy they cannot act as intermediate states.
A 2 Â 2 NOCI calculation in the basis of the two cluster basis states W S 0 S 1 and W S 1 S 0 gives the diabatic S [1] and S [2] states. The relative energies of these diabatic states are split by about 0.2 eV in the case of the pairs taken from stack A and stack B. The singlet excitation is delocalised (by symmetry) over the two molecules, as shown by the weights of the cluster basis states (see Table 3). In contrast, for the inter-stack pair the diabatic S[1] and S [2] states can be interpreted as being a localised singlet excitation on one of the molecules (see the weights in Table 3) and negligible energy splitting is obtained. The 1 TT state has a triplet excited state localised on each molecule, and in this model it is equal to the cluster basis state W1 TT .
To investigate the role of charge transfer states, a 4 Â 4 NOCI calculation in the basis of the cluster basis states W S 0 S 1 and W S 1 S 0 together with the two charge-transfer basis states W CT1 and W CT2 was performed. The relative energies and weights of these states are shown in Table 3. The charge-transfer basis states (weakly) interact with the diabatic S [1] and S [2] states, leading to a small energy lowering and to non-zero weights of the charge-transfer basis states. In contrast, the inclusion of the charge-transfer basis      Table 4). The magnitudes of the computed electronic couplings in both stacks A and B are sufficiently large that SF can occur efficiently [8]. The computed electronic couplings of the inter-stack pair are nearly close to zero. These results are not surprising since the band structure shows dispersion mainly along the stack direction of the crystal structure.

Conclusions
We have used a nonorthogonal configuration interaction approach for calculating the electronic coupling between the lowest diabatic excited singlet states and the 1 TT states. The diabatic S [1] and S [2] states can be interpreted as the first molecular singlet excitation, delocalised over two molecules. These states are indicative for delocalisation of the singlet excitation over the stack. These diabatic states (weakly) interact with the charge-transfer basis states. The diabatic 1 TT state can be interpreted as having a localised triplet excitation on each molecule. The computed electronic couplings between the diabatic S[1] and S [2] states and the 1 TT states are in the meV range, which is sufficient for the SF process. The inclusion of charge-transfer basis states enhances the computed electronic couplings and they act only as virtual states in the SF process. These results are obtained for antisymmetrized products of molecular CASSCF wavefunctions, which are found with a minimal active space consisting of only the two frontier molecular orbitals of the neutral molecule. The resulting couplings will of course change if more accurate molecular wavefunctions are employed. Nevertheless, the present results do indicate that this biradicaloid molecule is indeed a potential candidate as a SF chromophore. The nonorthogonal configuration interaction approach for calculating the electronic coupling between the diabatic excited singlet states and the 1 TT states is feasible and allows for a clear chemical interpretation of the diabatic states.