Communication: Active-space decomposition for molecular dimers

Interest in the excitonic properties of molecular crystals and aggregates is driven by their prominent role in established and emerging technologies such as organic light emitting devices^1,2 and cells,³ field effect transistors,⁴ photovoltaics,^5,6 and chemical sensors.⁷ Additionally, they are of biological interest for their role in photosynthesis,^8,9 and photodegradation and charge transport in DNA.^10–12 An accurate description of excitons in these systems is crucial to the further development of these materials and to the improved understanding of fundamental processes in biology.

The development of electronic structure methods for these excitonic systems has, however, proven challenging. Although density functional theory (DFT) is most commonly used, it has several well-known limitations including overly stabilized charge-transfer (CT) states, improper long-range exchange behavior, an undesirable dependence of electron transfer properties on the exchange-functional, and a lack of double excitations.¹³ Deficiencies in CT and long-range behaviors have been partially rectified by constrained DFT¹⁴ and range-separated hybrid functionals.¹⁵ The double excitations problem seems still open.¹⁶ The simplest wavefunction method for excited states, configuration interaction singles (CIS), overestimates the energy of CT states¹⁷ and does not describe double excitations. In equation-of-motion coupled cluster (EOM-CC) theory¹⁸ including its local variants,^19,20 triple displacements (the so-called EOM-CCSDT or CC3 models) are needed to quantitatively describe excited states with strong double excitation character, which are nonetheless expensive.

Electronic structure methods based on the active-space model have the potential to treat neutral, singly excited, doubly excited, and CT states on an equal footing. The complete-active-space method (CAS) and active-space methods with occupation restrictions^21–25 are well suited for calculations of excitonic properties, but their factorial scaling of computational cost with respect to number of active orbitals makes calculations of excitons in molecular aggregates computationally demanding (if even possible). We attempt to alleviate this high computational demand by tailoring an active-space model especially for molecular dimers. The dimer model is immediately applicable to processes such as charge or energy transfer and exciton fission. Furthermore, our model could be extended to treat general oligomers.

The theory developed in this work is based on tensor decomposition using molecular geometry. Tensor decomposition is a common concept in modern electronic structure theories. For instance, density fitting^26–28 and its extensions^29,30 approximate 2-electron integrals as a linear combination of products of low-rank integrals. Tensor factorization has been introduced into the so-called local electron correlation methods.³¹ Density matrix renormalization group (DMRG)^32,33 and tensor product state wavefunctions^34–37 are a compact representation of strongly correlated wavefunctions. Recently, low-rank representation in finite element wavefunctions have also been studied.³⁸ Our model combines the tensor decomposition concept with physical motivations from the wavefunctions of molecular aggregates.

In our theory a dimer wavefunction is parameterized by

where A and B label fragments (monomer units), I and J label fragment states (i.e., functions that diagonalize a fragment Hamiltonian) and K labels dimer states. The wavefunction is implicitly antisymmetrized with respect to the interchange of electrons. Although we focus on homodimers, the fragments A and B can be different molecules with different organizations of active spaces. This decomposition is exact when I and J run over a complete set of fragment states as we will demonstrate numerically. In practice, the number of states that need to be included in simulations of molecular dimers will be shown to be fairly small. The key result of this work is that using this ansatz, one can diagonalize the dimer Hamiltonian without forming the direct product configurations explicitly, leading to an algorithm that only requires intermediates that scale linearly with respect to the number of configurations in a monomer calculation (N_conf). Therefore, active spaces of the monomer in dimer calculations can be as large as that used in monomer calculations. This is in contrast to the quasi-complete-active-space (QCAS) method²² that uses a similar wavefunction ansatz. Our approach is similar in spirit to Molecular Orbital Valence Bond (MOVB)³⁹ and constrained DFT configuration interaction (CDFT-CI) methods in that we construct adiabatic states from a set of states that have some desired local diabatic character.^14,40 An important distinction with these methods is that we explicitly construct our states to be orthogonal.

First, we present an algorithm to compute Hamiltonian matrix elements between

where ζ is a combination of all the indices of orbitals belonging to subspace

in which ρ and σ denote spins, and p, q, and s label active orbitals of monomer A, and

Since the dimer basis functions in Eq. (1) and the Hamiltonian [Eq. (2)] are expressed in terms of direct products, the matrix elements can be computed by

using a monomer intermediate

We have applied this method to the stacked and T-shaped benzene dimer configurations as well as the slipped-parallel naphthalene dimer configuration. As a first test case, we consider a benzene dimer in a stacked geometry, as shown in Fig. 2(a). We start by examining the convergence of our method to a complete-active-space configuration interaction (CAS-CI) calculation in the full dimer space using the cc-pVDZ basis set.⁴⁴ Each monomer has an active space with 6 electrons and 6 orbitals so that the dimer has an active space of 12 electrons in 12 orbitals. The dimer active orbitals are determined automatically by taking those dimer orbitals that have maximum overlap with the active orbitals of isolated monomer fragments. Then the dimer active orbitals are localized using the Pipek–Mezey procedure⁴⁵ (by rotating among the 12 dimer active orbitals to maximize the square of the Mulliken populations of these orbitals on each atom). The overlap between these localized orbitals and the active orbitals of the isolated monomer fragments is used to assign each orbital to a particular fragment. Next, we compute diabatic monomer states using CAS-CI with the localized monomer active spaces. For example, to compute the “singlet” states of monomer A, a CAS-CI calculation is performed in which the active space contains monomer A’s localized active orbitals and the closed space includes the closed orbitals of monomer A and the filled orbitals of monomer B, i.e., a CAS-CI calculation of monomer A where monomer B is included as an embedding potential.

In our method, so far, one needs to specify what kinds of monomer states should be included and how many. In the simplest case, only monomer “singlet” (S) states are included. Our code can also compute matrix elements between arbitrary dimer states, such as CT states and different monomer spin states [“triplet” (T), “quintet” (Q), etc.]. The high-spin wavefunctions of monomer states (m_S = S) are computed via a standard full configuration interaction (FCI) algorithm⁴⁶ and the rest of the spin multiplet is generated by successive applications of the spin lowering operator,

Figure 3 demonstrates convergence of both the total ground state energy and the splitting between the first two excited singlet states for a stacked benzene dimer with d = 4.0  Å [see Fig. 2(a)] using different compositions and sizes of the monomer state bases. They are compared to a CAS-CI calculation in the full dimer space. As more states are included in the monomer subspaces, the ground state converges to the result obtained by CAS-CI in the full dimer active space. Interestingly, we note that in both cases CT-like monomer states are required to approach the benchmark result. This indicates the importance of polarization in the ground and excited states of the benzene stack. On the other hand, the calculations including T states are indistinguishable from those without.

We next consider the potential energy surfaces of the lowest-lying excited singlet electronic states of the benzene stack as a function of the separation distance, d, comparing the total energies from a full CAS-CI calculation to the total energies from our method using 16 states in the S, T, and CT monomer subspaces (Fig. 4). All surfaces computed with our method match the surfaces computed using a full CAS-CI approach to within 30 meV (≈1 mE_h) in the shaded region denoting the range of computed equilibrium distances^43,47 and even at least qualitatively at the extremely small separation of about 2.5 Å. This is especially encouraging considering that a typical calculation with our method required a total of 21 s (2.3 s for localization, 17.1 s computing localized CAS-CI states, 0.9 s forming the dimer Hamiltonian, and 0.2 s performing a Davidson diagonalization of the dimer Hamiltonian) on one processor compared to about 190 s for a comparable CAS-CI calculation in the full dimer space (9 states with 5 Davidson iterations). This speedup will become even more dramatic for larger active spaces due to the factorial scaling of the size of the CI vector with respect to the number of active orbitals.

We repeated the same analysis for the T-shaped benzene dimer. Since the spatial overlap between the two active spaces in the T-shaped dimer is smaller than that in the stacked dimer, the total and relative energies for the T-shaped dimer converged to the CAS-CI limit faster than they did for the stacked dimer. The numerical results are compiled in the supplementary material.⁴¹ 

In addition, we investigated the naphthalene dimer. With 20 electrons and 20 orbitals in the full π-space, the naphthalene dimer has 5.92 × 10⁹ configuration state functions and 3.41 × 10¹⁰ determinants, rendering conventional complete-active-space methods almost intractable; storing a single determinant-based CI vector, for example, would require over 250 GB. With our decomposition strategy, however, the full complete active space is never formed. We explore the convergence of the first two excited states of the naphthalene dimer in the slipped-parallel equilibrium configuration reported in Ref. 48 and shown in Fig. 2(b). The results for the first two excited states are shown in Fig. 5 using the cc-pVDZ basis set. Again, we see that CT states are important to the description of the excited states, while T states are less important. Calculations with S, CT, T subspaces and

In conclusion, we have presented here a physically motivated dimer wavefunction ansatz, in which the dimer wavefunction is represented in a basis of direct products of monomer wavefunctions. While this ansatz is quite general, we considered specifically the computation of complete-active-space wavefunctions for molecular dimers while only constructing monomer CAS-CI wavefunctions. This method will allow us to perform CAS-CI calculations on dimers that would be otherwise prohibitive.

Our method will be most powerful in situations where primarily diabatic states are of interest. These include electronic energy transfer processes where localized excited states are needed, charge transfer processes where charged diabatic states are needed, and singlet fission and triplet-triplet annihilation processes where (localized) multi-excitonic states are needed. This is especially true for molecular dimers that would be too large to be treated with a conventional CAS approach.

Possible extensions to our method are myriad. The approach could be modified to treat general oligomers by constructing the total wavefunction from direct products of N fragment wavefunctions (rather than 2). Orbital optimization could be performed by computing density matrices of the dimer space. Dynamic correlation could be introduced with a perturbative approach. Finally, we could improve the convergence of the method with respect to number of monomer states by using different criteria to choose the monomer states in which the dimer states are expanded. All of these possibilities are currently under investigation in our group.

T.Sh. has been supported by the Initiative for Sustainability and Energy at Northwestern (ISEN) Booster Award and the start-up fund provided by Northwestern University. T.Se. is grateful to the Department of Energy (DOE) (Grant No. DE-FG02-04ER15612) for support. S.M.P. is a DOE Office of Science Graduate Fellow.