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 μEh 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 shown3 that by choosing monomer wave functions that diagonalize a monomer Hamiltonian within the orthogonal subspace (i.e.,
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 CIJK⋯ 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 algorithm17,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
Comparison of wave functions in ASD-DMRG and in conventional ab initio DMRG.
The chain is initiated by computing a user-defined number of low-lying eigenstates of the first site. We will represent the ith single site as •i. These states are then renormalized to form a block, i.e., L1 ← •1 with Li referring to a block containing i sites. Next, a new site is added to the chain to form L1•2, and several low-lying eigenstates of this block-site system are computed and renormalized as L2 ← L1•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 Li−1•iRn−i is diagonalized, which contains a left block describing i − 1 sites Li−1, a right block describing n − i sites Rn−i, and the ith 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,
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 Li ← Li−1•i.
To diagonalize the Hamiltonian during the sweep, we resolve the Hamiltonian in terms of the renormalized block states
in which we have introduced the following operators:
Here X is either B or i,
where Z is either L or R, λ is either μ or ν, and
where
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 White19 to the RDM [Eq. (5)], in which the RDM is supplemented by
where a > 0 is an arbitrary constant, and
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 D6h symmetry using DFT with B3LYP functional23,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(6n,6n). 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 space26 indicate that ASD-DMRG with M = 128 recovered the CAS configuration interaction to within 10−7Eh. 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 × 1016 total configurations, is intractable even for today's most powerful supercomputers.
Total energies computed using M = 128 and energy differences obtained with selected smaller values of M for stacks of n benzene molecules (in Eh).
. | . | ΔE . | |$\smash{\raise -2pt\hbox{{E}_{\rm 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}}}$| . | ||||
---|---|---|---|---|---|---|---|
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 D2h symmetry with the def2-TZVPP basis set22 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 × 107 and 2 × 1012, 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 mEh. 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 mEh. The trimer calculation with M = 48 required ca. 1500 s per sweep using 256 CPU cores, although our program can be further optimized.
Total energies (in Eh) obtained for the PDI dimer and trimer at selected values of M.
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 mEh 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 adaptation28,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).