We extend the active space decomposition method, recently developed by us, to more than two active sites using the density matrix renormalization group algorithm. The fragment wave functions are described by complete or restricted active-space wave functions. Numerical results are shown on a benzene pentamer and a perylene diimide trimer. It is found that the truncation errors in our method decrease almost exponentially with respect to the number of renormalization states *M*, allowing for numerically exact calculations (to a few μ*E*_{h} or less) with *M* = 128 in both cases. This rapid convergence is because the renormalization steps are used only for the interfragment electron correlation.

The electronic structure of organic aggregates and crystals has been studied to characterize and predict electron and exciton dynamics in organic photovoltaics and light emitting devices.^{1} While the majority of computational studies of such processes has been performed using density functional theory (DFT), processes that involve multiple excitons, such as singlet fission,^{2} cannot be well described by standard time-dependent DFT, which has stimulated the development of accurate wave function theories for these systems.

We have recently developed a method, called active-space decomposition (ASD),^{3} which alleviates the large computational cost of active-space methods by tailoring the wave function ansatz to molecular dimers. In ASD, a dimer wave function Ψ is compactly expressed as a linear combination of tensor products of monomer wave functions, i.e.,

where *A* and *B* label monomers, and *I* and *J* label orthonormal monomer states. This expression is exact when the complete set of states is included in the sum; practically the number of states in the summation is set to a finite value. We have shown^{3} that by choosing monomer wave functions that diagonalize a monomer Hamiltonian within the orthogonal subspace (i.e.,

^{4}ASD is compatible with any wave function method for which transition density matrices are available. In particular, efficient implementation of ASD has been realized with complete active space (CAS) and restricted active space

^{5}(RAS) monomer wave functions. ASD has been used to compute model Hamiltonians for dynamical processes such as singlet fission in tetracene and pentacene dimers

^{6}and charge and triplet energy transfer in benzene dimers.

^{4}

Here, we extend ASD beyond dimers to incorporate an arbitrary number of fragments

where *I*, *J*, and *K* formally run over the complete set of states on each fragment. The coefficient tensor *C*_{IJK⋯} is approximated by matrix product states used in the density matrix renormalization group (DMRG) algorithm,

in which the dimension of matrices **C** is set to a predetermined value *M*. First described by White and used to compute ground and excited states of the Heisenberg spin chain,^{7,8} DMRG has proven to be a powerful tool that provides numerically exact solutions with polynomial cost even for strongly correlated systems.^{9–14} In chemical applications, DMRG has enabled highly accurate calculations on strongly correlated electronic structure of molecules or complexes that are intractable with traditional active-space methods.^{15,16}

We start by sketching the ASD-DMRG algorithm, which is formulated in terms of a one-dimensional chain of *n* sites (see Fig. 1). In the conventional *ab initio* DMRG, these sites are chosen to be one-electron orbitals.^{9} In this work, we choose each site to be the CAS or RAS wave function of a single molecule or fragment, alleviating the orbital ordering problem in the conventional *ab initio* DMRG algorithm^{17,18} and allowing for applications to three-dimensional structures. As a consequence, the dimensionality of the single site, *d*, is much larger in ASD-DMRG than in conventional approaches. However, a key result of this study is that this allows *M* to be very small. Since calculations involving the full configurational space of a molecular dimer are intractable for large molecules of materials interest, we have implemented the single-site variant of the DMRG algorithm.^{19} The single-site algorithm converges more smoothly than the two-site algorithm, though it is more prone to getting stuck in local minima during the optimization (see below).^{19,20}

The chain is initiated by computing a user-defined number of low-lying eigenstates of the first site. We will represent the *i*th single site as •_{i}. These states are then renormalized to form a block, i.e., *L*_{1} ← •_{1} with *L*_{i} referring to a block containing *i* sites. Next, a new site is added to the chain to form *L*_{1}•_{2}, and several low-lying eigenstates of this block-site system are computed and renormalized as *L*_{2} ← *L*_{1}•_{2}. This process is continued until each site has been added to the chain. During this growing phase, the sites that have yet to be incorporated into the block are included at the mean-field level. This is equivalent to using a tensor product of the Hartree–Fock ground state of each site as the initial guess. Such a mean-field initial guess is well-defined in the case where each site is a distinct molecule, although other initial guesses may be necessary in the case where sites are covalently bound.

The final stage of the algorithm involves “sweeping” along the chain to optimize the renormalized states. In each step of the sweep, a system *L*_{i−1}•_{i}*R*_{n−i} is diagonalized, which contains a left block describing *i* − 1 sites *L*_{i−1}, a right block describing *n* − *i* sites *R*_{n−i}, and the *i*th single site, •_{i}. The wave function at this step is written as

where |ϕ^{i}⟩ is a Slater determinant only with orbitals on the site *i*,

*L*

_{i−1}•

_{i}part of the chain can be obtained by the Schmidt decomposition of the reduced density matrix (RDM),

where the trace is over all states in the right block. Thus, one computes and diagonalizes the RDM in each step of the sweep and uses the *M* eigenvectors corresponding to the largest eigenvalues to renormalize *L*_{i} ← *L*_{i−1}•_{i}.

To diagonalize the Hamiltonian during the sweep, we resolve the Hamiltonian in terms of the renormalized block states

^{11}we write the Hamiltonian as

in which we have introduced the following operators:

Here *X* is either *B* or *i*,

*p*on site/block

*X*,

*h*

_{pq}is a matrix element of the one-electron operator, and (

*pq*|

*rs*) is a matrix element of the two-electron Coulomb operator. In our algorithm, we compute and store transition density matrices whose indices belong to the same block using an algorithm developed in Refs. 3 and 4. For instance,

where *Z* is either *L* or *R*, λ is either μ or ν, and

where

*L*

_{μ}.

Active orbitals are identified as follows. We start by picking orbitals from monomer Hartree–Fock calculations using minimal basis. Then we perform Hartree–Fock on the full system and localize the canonical molecular orbitals to fragments using the Pipek–Mezey algorithm.^{21} Finally, overlaps between these fragment-localized orbitals and the minimal-basis orbitals are used to automatically select active orbitals on each site.

As is commonly known,^{19,20} the one-site DMRG algorithm is susceptible to getting stuck in local minima of the parameter space due to the fact that the algorithm does not spontaneously expand the range of quantum numbers in each site or block. For example, if a poor initial guess for the chain includes only neutral fragments and the total charge is constrained to be neutral, then the algorithm will keep only neutral fragment states although charge transfer configurations may be important in the exact ground state. To remedy this problem, we use the perturbative correction of White^{19} to the RDM [Eq. (5)], in which the RDM is supplemented by

where *a* > 0 is an arbitrary constant, and

*a*is used to control the strength of the perturbation during the sweep. We typically start with a value of 10

^{−3}and decrease it by one order of magnitude when the energy from successive sweeps changes by less than 10

^{−7}

*E*

_{h}. When

*a*reaches a minimum value (

*a*< 10

^{−7}in our calculations) the perturbation is turned off. With this scheme, we have seen that the final converged energies have no dependence on the initial guess.

Next, we show the numerical results on the ground state of a benzene stack consisting of 2–5 molecules with the def2-SVP basis set.^{22} The molecules are arranged cofacially with a separation of 4.0 Å (see Fig. 2(a)). The geometry of each benzene molecule was obtained in *D*_{6h} symmetry using DFT with B3LYP functional^{23,24} and the def2-TZVPP basis set.^{22} Geometry optimization was performed using TURBOMOLE.^{25} Each benzene molecule is a single site containing six π orbitals and six electrons. Using CAS wave functions on each site, the full stack is the equivalent of CAS(6*n*,6*n*). Calculations on a benzene dimer [for which the configuration interaction in CAS(12,12) is tractable] and those on a benzene trimer with a smaller active space^{26} indicate that ASD-DMRG with *M* = 128 recovered the CAS configuration interaction to within 10^{−7}*E*_{h}. Thus, we consider *M* = 128 to be fully converged for all benzene stack calculations. In Table I, the total energies are compared to those obtained using smaller values of *M*. It is worth noting that the errors decrease almost exponentially with respect to *M*, and with as small as *M* = 32, numerically exact results were obtained. The pentamer calculation with *M* = 128 required about 2200 s per sweep on 128 CPU cores (9 sweeps), whereas the corresponding CAS(30,30) calculation, with 2.4 × 10^{16} total configurations, is intractable even for today's most powerful supercomputers.

. | . | ΔE
. | |$\smash{\raise -2pt\hbox{{E}_{\rm tot}}}$|$E tot $ . | ||||
---|---|---|---|---|---|---|---|

n
. | Active space . | M = 1
. | M = 4
. | M = 8
. | M = 16
. | M = 32
. | M = 128
. |

2 | CAS(12,12) | 0.000501 | 0.000168 | 0.000134 | 0.000069 | < 0.000001 | −461.179327 |

3 | CAS(18,18) | 0.001018 | 0.000343 | 0.000274 | 0.000141 | 0.000001 | −691.768264 |

4 | CAS(24,24) | 0.001537 | 0.000520 | 0.000415 | 0.000214 | 0.000002 | −922.357173 |

5 | CAS(30,30) | 0.002055 | 0.000696 | 0.000557 | 0.000287 | 0.000002 | −1152.946072 |

. | . | ΔE
. | |$\smash{\raise -2pt\hbox{{E}_{\rm tot}}}$|$E tot $ . | ||||
---|---|---|---|---|---|---|---|

n
. | Active space . | M = 1
. | M = 4
. | M = 8
. | M = 16
. | M = 32
. | M = 128
. |

2 | CAS(12,12) | 0.000501 | 0.000168 | 0.000134 | 0.000069 | < 0.000001 | −461.179327 |

3 | CAS(18,18) | 0.001018 | 0.000343 | 0.000274 | 0.000141 | 0.000001 | −691.768264 |

4 | CAS(24,24) | 0.001537 | 0.000520 | 0.000415 | 0.000214 | 0.000002 | −922.357173 |

5 | CAS(30,30) | 0.002055 | 0.000696 | 0.000557 | 0.000287 | 0.000002 | −1152.946072 |

One of the advantages of ASD-DMRG over conventional DMRG is that occupation restrictions can be introduced on individual sites. To show this, we computed the ground state of the perylene diimide (PDI) dimer and trimer using RAS wave functions for each molecule (def2-SVP). PDI crystallizes into a slip-stacked configuration as shown in Fig. 2(b).^{27} For monomers, we used geometries optimized using DFT under *D*_{2h} symmetry with the def2-TZVPP basis set^{22} and the B3LYP density functional.^{23,24} The dimer and trimer geometries were generated by adding a constant offset to the base unit of (−3.34, 1.11, 3.38) Å where *x* is the long axis of the molecule, *y* is the short axis of the molecule, and *z* is normal to the plane of the molecule. For single-site wave functions, we use RAS with 7 orbitals in RAS I, 4 orbitals in RAS II, and 7 orbitals in RAS III and a maximum of one hole in RAS I and one particle in RAS III. The tensor product of multiple RAS wave functions corresponds to a more general set of occupation restrictions than RAS itself: Each site can have one hole (particle) in RAS I (RAS III) rather than there being one hole (particle) allowed for the whole stack. The dimensionality of the product wave functions for the dimer and trimer is about 1 × 10^{7} and 2 × 10^{12}, respectively.

The total energies computed for the PDI dimer and trimer are shown in Table II. For the PDI dimer, we used up to *M* = 128. As in the benzene case, the convergence of the total energy with respect to *M* appears to be exponential. On the basis of these calculations we chose *M* = 48 to compute the ground state of the PDI trimer. The energy difference between *M* = 128 and *M* = 48 for the dimer is about 0.06 m*E*_{h}. Noting that error in the total energy of the benzene stacks for a fixed *M* was roughly proportional to *n* − 1, we expect the total energy for the PDI trimer at *M* = 48 to be accurate to within 0.2 m*E*_{h}. The trimer calculation with *M* = 48 required ca. 1500 s per sweep using 256 CPU cores, although our program can be further optimized.

M
. | Dimer . | Trimer . |
---|---|---|

1 | −2644.328613 | −3966.487183 |

4 | −2644.331085 | −3966.492073 |

8 | −2644.331681 | −3966.493276 |

16 | −2644.331990 | −3966.493900 |

32 | −2644.332199 | −3966.494332 |

48 | −2644.332283 | −3966.494512 |

72 | −2644.332327 | … |

96 | −2644.332339 | … |

128 | −2644.332343 | … |

M
. | Dimer . | Trimer . |
---|---|---|

1 | −2644.328613 | −3966.487183 |

4 | −2644.331085 | −3966.492073 |

8 | −2644.331681 | −3966.493276 |

16 | −2644.331990 | −3966.493900 |

32 | −2644.332199 | −3966.494332 |

48 | −2644.332283 | −3966.494512 |

72 | −2644.332327 | … |

96 | −2644.332339 | … |

128 | −2644.332343 | … |

In conclusion, we have generalized the ASD method to incorporate multiple fragments (i.e., active sites). In this work, the electronic structure of each fragment has been described by CAS and RAS wave functions. Our method is based on the matrix-product approximation of the coefficients [Eq. (3)], in which the wave functions are optimized using the density matrix renormalization group algorithm. We have shown that results with chemical accuracy (i.e., 1 m*E*_{h} for total energies) can be obtained with a fairly small number of renormalization states *M*. This rapid convergence of the ASD-DMRG energy with respect to *M* is the result of breaking the system into physically motivated fragments such that the renormalization steps are used to describe the interchromophore electron correlation, and not the intrachromophore correlation. Spin adaptation^{28,29} has not been done in this work. Application of ASD-DMRG to covalently bonded systems and investigation of the effects of site ordering will be performed in the future. All the programs reported in this work have been implemented in the open-source BAGEL package.^{30}

This work has been supported by Office of Basic Energy Sciences, U.S. Department of Energy (Grant No. DE-FG02-13ER16398).