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.,

(1)

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.,

$\hat{H}^A|\Phi ^A_I\rangle = E^A_I|\Phi ^A_I\rangle$
ĤA|ΦIA=EIA|ΦIA⁠), dimer wave functions rapidly converge to the exact solution with respect to the number and type of monomer states. This product basis is also related to quasi-diabatic representation of the model space.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 space5 (RAS) monomer wave functions. ASD has been used to compute model Hamiltonians for dynamical processes such as singlet fission in tetracene and pentacene dimers6 and charge and triplet energy transfer in benzene dimers.4 

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

(2)

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,

(3)

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

FIG. 1.

Comparison of wave functions in ASD-DMRG and in conventional ab initio DMRG.

FIG. 1.

Comparison of wave functions in ASD-DMRG and in conventional ab initio DMRG.

Close modal

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 L12, and several low-lying eigenstates of this block-site system are computed and renormalized as L2L12. 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−1iRni is diagonalized, which contains a left block describing i − 1 sites Li−1, a right block describing ni sites Rni, and the ith single site, •i. The wave function at this step is written as

(4)

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

$|\Phi ^{L_{i-1}}_\mu \rangle$
|ΦμLi1 is a renormalized left-block state,
$|\Phi ^{R_{n-i}}_\nu \rangle$
|ΦνRni
is a renormalized right-block state, and
$|\Phi ^i_{\mu \nu }\rangle \equiv \sum _\phi C^{i,\phi }_{\mu \nu } |\phi^i\rangle$
|ΦμνiϕCμνi,ϕ|ϕi
. The key insight behind the DMRG algorithm is that the optimal block states to represent the combined Li−1i part of the chain can be obtained by the Schmidt decomposition of the reduced density matrix (RDM),

(5)

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 LiLi−1i.

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

$\lbrace \Phi ^L_\mu \rbrace$
{ΦμL} and
$\lbrace \Phi ^R_\nu \rbrace$
{ΦνR}
. Defining
$|\Phi ^B_{\mu \nu }\rangle \equiv |\Phi ^L_\mu \rangle \otimes |\Phi ^R_\nu \rangle$
|ΦμνB|ΦμL|ΦνR
and following the notation of Kurashige and Yanai,11 we write the Hamiltonian as

(6)

in which we have introduced the following operators:

(7)
(8)
(9)
(10)
(11)

Here X is either B or i,

$\hat{p}_X$
p̂X (
$\hat{p}_X^{\dagger })$
p̂X)
is the annihilation (creation) operator for orbital p on site/block X, hpq 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,

(12)

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

$\Gamma ^{Z,\lambda ^{\prime }\lambda }_{p^\dagger q^\dagger r}$
ΓpqrZ,λλ⁠,
$\Gamma ^{Z,\lambda ^{\prime }\lambda }_{p q}$
ΓpqZ,λλ
, and
$\Gamma ^{Z,\lambda ^{\prime }\lambda }_{p}$
ΓpZ,λλ
are likewise defined. The block operators of Eq. (6) are resolved using these transition densities, e.g.,

(13)

where

$N_{L_\mu }$
NLμ is the number of electrons in 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 White19 to the RDM [Eq. (5)], in which the RDM is supplemented by

(14)

where a > 0 is an arbitrary constant, and

$\hat{O}^i_\tau$
Ôτi is a site Hamiltonian operator (we use
$\hat{p}_i$
p̂i
,
$\hat{p}_i^\dagger$
p̂i
, and
$\hat{p}_i^\dagger \hat{q}_i$
p̂iq̂i
). The constant 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−7Eh. 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 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.

FIG. 2.

Geometries of the (a) benzene and (b) PDI stacks used in this work.

FIG. 2.

Geometries of the (a) benzene and (b) PDI stacks used in this work.

Close modal
Table I.

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}}}$|E tot
nActive spaceM = 1M = 4M = 8M = 16M = 32M = 128
CAS(12,12) 0.000501 0.000168 0.000134 0.000069 < 0.000001 −461.179327 
CAS(18,18) 0.001018 0.000343 0.000274 0.000141 0.000001 −691.768264 
CAS(24,24) 0.001537 0.000520 0.000415 0.000214 0.000002 −922.357173 
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
nActive spaceM = 1M = 4M = 8M = 16M = 32M = 128
CAS(12,12) 0.000501 0.000168 0.000134 0.000069 < 0.000001 −461.179327 
CAS(18,18) 0.001018 0.000343 0.000274 0.000141 0.000001 −691.768264 
CAS(24,24) 0.001537 0.000520 0.000415 0.000214 0.000002 −922.357173 
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.

Table II.

Total energies (in Eh) obtained for the PDI dimer and trimer at selected values of M.

MDimerTrimer
−2644.328613 −3966.487183 
−2644.331085 −3966.492073 
−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 … 
MDimerTrimer
−2644.328613 −3966.487183 
−2644.331085 −3966.492073 
−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).

1.
V.
Coropceanu
,
J.
Cornil
,
D. A.
da Silva Filho
,
Y.
Olivier
,
R.
Silbey
, and
J.-L.
Brédas
,
Chem. Rev.
107
,
926
(
2007
).
2.
M. B.
Smith
and
J.
Michl
,
Chem. Rev.
110
,
6891
(
2010
).
3.
S. M.
Parker
,
T.
Seideman
,
M. A.
Ratner
, and
T.
Shiozaki
,
J. Chem. Phys.
139
,
021108
(
2013
).
4.
S. M.
Parker
and
T.
Shiozaki
,
J. Chem. Theory Comput.
10
,
3738
(
2014
).
5.
J.
Olsen
,
B. O.
Roos
,
P.
Jørgensen
, and
H. J. Aa.
Jensen
,
J. Chem. Phys.
89
,
2185
(
1988
).
6.
S. M.
Parker
,
T.
Seideman
,
M. A.
Ratner
, and
T.
Shiozaki
,
J. Phys. Chem. C
118
,
12700
(
2014
).
7.
8.
S. R.
White
,
Phys. Rev. B
48
,
10345
(
1993
).
9.
S. R.
White
and
R. L.
Martin
,
J. Chem. Phys.
110
,
4127
(
1999
).
10.
D.
Zgid
and
G. K.-L.
Chan
,
Annu. Rep. Comput. Chem.
5
,
149
(
2009
).
11.
Y.
Kurashige
and
T.
Yanai
,
J. Chem. Phys.
130
,
234114
(
2009
).
12.
K. H.
Marti
and
M.
Reiher
,
Z. Phys. Chem.
224
,
583
(
2010
).
13.
G. K.-L.
Chan
and
S.
Sharma
,
Annu. Rev. Phys. Chem.
62
,
465
(
2011
).
14.
S.
Wouters
and
D.
Van Neck
,
Eur. Phys. J. D
68
,
272
(
2014
).
15.
Y.
Kurashige
,
G. K.-L.
Chan
, and
T.
Yanai
,
Nature Chem.
5
,
660
(
2013
).
16.
S.
Sharma
,
K.
Sivalingam
,
F.
Neese
, and
G. K.-L.
Chan
,
Nature Chem.
6
,
927
(
2014
).
17.
G.
Moritz
,
B. A.
Hess
, and
M.
Reiher
,
J. Chem. Phys.
122
,
024107
(
2005
).
18.
J.
Rissler
,
R. M.
Roack
, and
S. R.
White
,
Chem. Phys.
323
,
519
(
2006
).
19.
S. R.
White
,
Phys. Rev. B
72
,
180403
(
2005
).
20.
D.
Zgid
and
M.
Nooijen
,
J. Chem. Phys.
128
,
144116
(
2008
).
21.
J.
Pipek
and
P. G.
Mezey
,
J. Chem. Phys.
90
,
4916
(
1989
).
22.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
23.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
24.
C.
Lee
,
W.
Yang
, and
R.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
25.
F.
Furche
,
R.
Ahlrichs
,
C.
Hättig
,
W.
Klopper
,
M.
Sierka
, and
F.
Weigend
,
WIREs Comput. Mol. Sci.
4
,
91
(
2014
).
26.
See supplementary material at http://dx.doi.org/10.1063/1.4902991 for numerical tests of our implementation on the benzene dimer and trimer using a CAS(4,4) active space.
27.
O.
Guillermet
,
M.
Mossoyan-Déneux
,
M.
Giorgi
,
A.
Glachant
, and
J. C.
Mossoyan
,
Thin Solid Films
514
,
25
(
2006
).
28.
D.
Zgid
and
M.
Nooijen
,
J. Chem. Phys.
128
,
014107
(
2008
).
29.
S.
Sharma
and
G. K.-L.
Chan
,
J. Chem. Phys.
136
,
144105
(
2012
).
30.
BAGEL, Brilliantly Advanced General Electronic-structure Library, See http://www.nubakery.org under the GNU General Public License.

Supplementary Material