Skip to Main Content

Skip Nav Destination

Article navigation

Rapid Communication| December 03 2014

Communication: Active space decomposition with multiple sites: Density matrix renormalization group algorithm

Shane M. Parker

0000-0002-1110-3393

;

Shane M. Parker

Department of Chemistry,

Northwestern University

, 2145 Sheridan Rd., Evanston, Illinois 60208,

USA

Search for other works by this author on:

PubMed

Google Scholar

Toru Shiozaki

Toru Shiozaki

Department of Chemistry,

Northwestern University

, 2145 Sheridan Rd., Evanston, Illinois 60208,

USA

Search for other works by this author on:

PubMed

Google Scholar

Author & Article Information

Shane M. Parker

^{^{https://orcid.org/0000-0002-1110-3393}}

,

Toru Shiozaki

Department of Chemistry,

Northwestern University

, 2145 Sheridan Rd., Evanston, Illinois 60208,

USA

J. Chem. Phys. 141, 211102 (2014)

https://doi.org/10.1063/1.4902991

CHORUS

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.¹ 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,² 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),³ 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.,

\begin{equation}|\Psi \rangle = {\sum _{IJ}}C_{IJ}\big|\Phi ^A_I\big\rangle \otimes \big|\Phi ^B_J\big\rangle \equiv \sum _{IJ} C_{IJ} \big|\Phi ^A_I\Phi ^B_J\big\rangle ,\end{equation}

| Ψ ⟩ = \sum_{I J} C_{I J} | Φ_{I}^{A} ⟩ \otimes | Φ_{J}^{B} ⟩ \equiv \sum_{I J} C_{I J} | Φ_{I}^{A} Φ_{J}^{B} ⟩,

(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 shown³ 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$

{\hat{H}}^{A} | Φ_{I}^{A} ⟩ = E_{I}^{A} | Φ_{I}^{A} ⟩

⁠), 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.⁴ 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⁵ (RAS) monomer wave functions. ASD has been used to compute model Hamiltonians for dynamical processes such as singlet fission in tetracene and pentacene dimers⁶ and charge and triplet energy transfer in benzene dimers.⁴

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

\begin{equation}|\Psi \rangle = \sum _{IJK\cdots }C_{IJK\cdots }\big|\Phi ^A_I\big\rangle \otimes \big|\Phi ^B_J\big\rangle \otimes \big|\Phi ^C_K\big\rangle \otimes \cdots ,\end{equation}

| Ψ ⟩ = \sum_{I J K \dots} C_{I J K \dots} | Φ_{I}^{A} ⟩ \otimes | Φ_{J}^{B} ⟩ \otimes | Φ_{K}^{C} ⟩ \otimes \dots,

(2)

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,

\begin{equation}C_{IJK\cdots } = \mathrm{Tr}[\mathbf {C}^{A,I}\mathbf {C}^{B,J}\mathbf {C}^{C,K}\cdots ],\end{equation}

C_{I J K \dots} = Tr [C^{A, I} C^{B, J} C^{C, K} \dots],

(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.⁹ 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.¹⁹ 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.

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

View large Download slide

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., L₁ ← •₁ with L_i referring to a block containing i sites. Next, a new site is added to the chain to form L₁•₂, and several low-lying eigenstates of this block-site system are computed and renormalized as L₂ ← L₁•₂. 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•_iR_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 ith single site, •_i. The wave function at this step is written as

\begin{align}|\Psi \rangle &= \sum ^M_{\mu \nu } \sum _\phi C^{i,\phi }_{\mu \nu }\big|\Phi ^{L_{i-1}}_\mu \big\rangle \otimes |\phi ^i\rangle \otimes \big|\Phi ^{R_{n-i}}_\nu \big\rangle \nonumber \\&= \sum ^M_{\mu \nu }\big|\Phi ^{L_{i-1}}_\mu \Phi ^i_{\mu \nu } \Phi ^{R_{n-i}}_\nu \big\rangle ,\end{align}

\begin{matrix} | Ψ ⟩ & = \sum_{μ ν}^{M} \sum_{ϕ} C_{μ ν}^{i, ϕ} | Φ_{μ}^{L_{i - 1}} ⟩ \otimes | ϕ^{i} ⟩ \otimes | Φ_{ν}^{R_{n - i}} ⟩ \\ = \sum_{μ ν}^{M} | Φ_{μ}^{L_{i - 1}} Φ_{μ ν}^{i} Φ_{ν}^{R_{n - i}} ⟩, \end{matrix}

(4)

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

$|\Phi ^{L_{i-1}}_\mu \rangle$

| Φ_{μ}^{L_{i - 1}} ⟩

is a renormalized left-block state,

$|\Phi ^{R_{n-i}}_\nu \rangle$

| Φ_{ν}^{R_{n - i}} ⟩

is a renormalized right-block state, and

$|\Phi ^i_{\mu \nu }\rangle \equiv \sum _\phi C^{i,\phi }_{\mu \nu } |\phi^i\rangle$

| Φ_{μ ν}^{i} ⟩ \equiv \sum_{ϕ} C_{μ ν}^{i, ϕ} | ϕ^{i} ⟩

⁠. The key insight behind the DMRG algorithm is that the optimal block states to represent the combined L_i−1•_i part of the chain can be obtained by the Schmidt decomposition of the reduced density matrix (RDM),

\begin{equation}\hat{\rho }_{L_{i-1}\bullet _i} = \text{Tr}\left[ |\Psi \rangle \langle \Psi |\right]_{R_{n-i}},\end{equation}

{\hat{ρ}}_{L_{i - 1} •_{i}} = Tr {[| Ψ ⟩ ⟨ Ψ |]}_{R_{n - i}},

(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 L_i ← L_i−1•_i.

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} ⟩ \equiv | Φ_{μ}^{L} ⟩ \otimes | Φ_{ν}^{R} ⟩

and following the notation of Kurashige and Yanai,¹¹ we write the Hamiltonian as

\begin{eqnarray}\hat{H}&= &\hat{H}^B + \hat{H}^i + \sum _{p_B}\big(\hat{T}^i_{p_B}\hat{p}_B + \hat{p}_B^\dagger \hat{T}^{i\dagger }_{p_B} \big) \nonumber \\&&+\, \sum _{p_i}\big(\hat{S}^B_{p_i}\hat{p}_i + \hat{p}_i^\dagger \hat{S}^{B\dagger }_{p_i} \big) + \sum _{p_iq_i}\hat{p}_i^\dagger \hat{q}_i\hat{Q}^B_{p_iq_i} \nonumber \\&&+ \, \sum _{p_ir_i}\big(\hat{p}_i^\dagger \hat{r}_i^\dagger \hat{P}^B_{p_ir_i} + \hat{P}^{B\dagger }_{p_ir_i}\hat{p}_i\hat{r}_i \big),\end{eqnarray}

\begin{matrix} \hat{H} & = & {\hat{H}}^{B} + {\hat{H}}^{i} + \sum_{p_{B}} ({\hat{T}}_{p_{B}}^{i} {\hat{p}}_{B} + {\hat{p}}_{B}^{†} {\hat{T}}_{p_{B}}^{i †}) \\ + \sum_{p_{i}} ({\hat{S}}_{p_{i}}^{B} {\hat{p}}_{i} + {\hat{p}}_{i}^{†} {\hat{S}}_{p_{i}}^{B †}) + \sum_{p_{i} q_{i}} {\hat{p}}_{i}^{†} {\hat{q}}_{i} {\hat{Q}}_{p_{i} q_{i}}^{B} \\ + \sum_{p_{i} r_{i}} ({\hat{p}}_{i}^{†} {\hat{r}}_{i}^{†} {\hat{P}}_{p_{i} r_{i}}^{B} + {\hat{P}}_{p_{i} r_{i}}^{B †} {\hat{p}}_{i} {\hat{r}}_{i}), \end{matrix}

(6)

in which we have introduced the following operators:

\begin{align}\hat{H}^X =& \sum _{p_Xq_X}\hat{p}_X^\dagger \hat{q}_X h_{p_Xq_X} + \frac{1}{2} \sum _{p_Xq_Xr_Xs_X}\hat{p}_X^\dagger \hat{q}_X^\dagger \hat{r}_X \hat{s}_X (p_Xs_X|q_Xr_X),\end{align}

\begin{matrix} {\hat{H}}^{X} = & \sum_{p_{X} q_{X}} {\hat{p}}_{X}^{†} {\hat{q}}_{X} h_{p_{X} q_{X}} + \frac{1}{2} \sum_{p_{X} q_{X} r_{X} s_{X}} {\hat{p}}_{X}^{†} {\hat{q}}_{X}^{†} {\hat{r}}_{X} {\hat{s}}_{X} (p_{X} s_{X} | q_{X} r_{X}), \end{matrix}

(7)

\begin{align}\hat{T}^i_{p_B} =& \sum _{q_ir_is_i}\hat{s}_i^\dagger \hat{r}_i^\dagger \hat{q}_i (p_Bs_i|q_ir_i),\end{align}

\begin{matrix} {\hat{T}}_{p_{B}}^{i} = & \sum_{q_{i} r_{i} s_{i}} {\hat{s}}_{i}^{†} {\hat{r}}_{i}^{†} {\hat{q}}_{i} (p_{B} s_{i} | q_{i} r_{i}), \end{matrix}

(8)

\begin{align}\hat{S}^B_{p_i} =& \sum _{q_B} \hat{q}_B^\dagger h_{p_iq_B} + \sum _{q_B,r_B,s_B}\hat{s}_B^\dagger \hat{r}_B^\dagger \hat{q}_B (p_is_B|q_Br_X),\end{align}

\begin{matrix} {\hat{S}}_{p_{i}}^{B} = & \sum_{q_{B}} {\hat{q}}_{B}^{†} h_{p_{i} q_{B}} + \sum_{q_{B}, r_{B}, s_{B}} {\hat{s}}_{B}^{†} {\hat{r}}_{B}^{†} {\hat{q}}_{B} (p_{i} s_{B} | q_{B} r_{X}), \end{matrix}

(9)

\begin{align}\hat{P}^B_{p_ir_i} =& \sum _{q_is_i}\hat{q}_i\hat{s}_i (p_Bs_i|q_ir_B),\end{align}

\begin{matrix} {\hat{P}}_{p_{i} r_{i}}^{B} = & \sum_{q_{i} s_{i}} {\hat{q}}_{i} {\hat{s}}_{i} (p_{B} s_{i} | q_{i} r_{B}), \end{matrix}

(10)

\begin{align}\hat{Q}^B_{p_iq_i} =& \sum _{r_is_i}\hat{r}_i^\dagger \hat{s}_i [ (p_Bq_B|r_is_i) - (p_Br_i|q_Bs_i) ].\end{align}

\begin{matrix} {\hat{Q}}_{p_{i} q_{i}}^{B} = & \sum_{r_{i} s_{i}} {\hat{r}}_{i}^{†} {\hat{s}}_{i} [(p_{B} q_{B} | r_{i} s_{i}) - (p_{B} r_{i} | q_{B} s_{i})] . \end{matrix}

(11)

Here X is either B or i,

$\hat{p}_X$

{\hat{p}}_{X}

(⁠

$\hat{p}_X^{\dagger })$

{\hat{p}}_{X}^{†})

is the annihilation (creation) operator for orbital 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,

\begin{align}\Gamma ^{Z,\lambda ^{\prime }\lambda }_{p^\dagger q} \equiv \big\langle \Phi ^Z_{\lambda ^{\prime }}\big|\hat{p}^\dagger \hat{q}\big|\Phi ^Z_\lambda \big\rangle ,\end{align}

\begin{matrix} Γ_{p^{†} q}^{Z, λ^{'} λ} \equiv ⟨Φ_{λ^{'}}^{Z} | {\hat{p}}^{†} \hat{q} | Φ_{λ}^{Z}⟩, \end{matrix}

(12)

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

$\Gamma ^{Z,\lambda ^{\prime }\lambda }_{p^\dagger q^\dagger r}$

Γ_{p^{†} q^{†} r}^{Z, λ^{'} λ}

⁠,

$\Gamma ^{Z,\lambda ^{\prime }\lambda }_{p q}$

Γ_{p q}^{Z, λ^{'} λ}

⁠, and

$\Gamma ^{Z,\lambda ^{\prime }\lambda }_{p}$

Γ_{p}^{Z, λ^{'} λ}

are likewise defined. The block operators of Eq. are resolved using these transition densities, e.g.,

\begin{eqnarray}&&\big\langle \Phi ^B_{\mu ^{\prime }\nu ^{\prime }}\big|\hat{P}^B_{p_ir_i}\big|\Phi ^B_{\mu \nu }\big\rangle \nonumber\\&&\quad = \delta _{\nu ^{\prime }\nu }\sum _{q_Ls_L}\Gamma ^{L,\mu ^{\prime }\mu }_{q_Ls_L} (p_is_L|r_iq_L) \nonumber \\&&\qquad +\, \delta _{\mu ^{\prime }\mu }\sum _{q_Rs_R}\Gamma ^{R,\nu ^{\prime }\nu }_{q_Rs_R} (p_is_R|r_iq_R) \nonumber \\&&\qquad -\, (-1)^{N_{L_\mu }}\sum _{q_R}\Gamma ^{R,\nu ^{\prime }\nu }_{q_R}\sum _{s_L}\Gamma ^{L,\mu ^{\prime }\mu }_{s_L} (p_is_L|r_iq_R) \nonumber \\&&\qquad +\, (-1)^{N_{L_\mu }}\sum _{q_L}\Gamma ^{L,\mu ^{\prime }\mu }_{q_L}\sum _{s_R}\Gamma ^{R,\nu ^{\prime }\nu }_{s_R} (p_is_R|r_iq_L),\qquad\end{eqnarray}

\begin{matrix} ⟨Φ_{μ^{'} ν^{'}}^{B} | {\hat{P}}_{p_{i} r_{i}}^{B} | Φ_{μ ν}^{B}⟩ \\ = δ_{ν^{'} ν} \sum_{q_{L} s_{L}} Γ_{q_{L} s_{L}}^{L, μ^{'} μ} (p_{i} s_{L} | r_{i} q_{L}) \\ + δ_{μ^{'} μ} \sum_{q_{R} s_{R}} Γ_{q_{R} s_{R}}^{R, ν^{'} ν} (p_{i} s_{R} | r_{i} q_{R}) \\ - {(- 1)}^{N_{L_{μ}}} \sum_{q_{R}} Γ_{q_{R}}^{R, ν^{'} ν} \sum_{s_{L}} Γ_{s_{L}}^{L, μ^{'} μ} (p_{i} s_{L} | r_{i} q_{R}) \\ + {(- 1)}^{N_{L_{μ}}} \sum_{q_{L}} Γ_{q_{L}}^{L, μ^{'} μ} \sum_{s_{R}} Γ_{s_{R}}^{R, ν^{'} ν} (p_{i} s_{R} | r_{i} q_{L}), \end{matrix}

(13)

where

$N_{L_\mu }$

N_{L_{μ}}

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.²¹ 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¹⁹ to the RDM [Eq. (5)], in which the RDM is supplemented by

\begin{align}\hat{\rho }_{L_{i-1}\bullet _i}^{\prime } = \hat{\rho }_{L_{i-1}\bullet _i} + a \sum _\tau \hat{O}^{i}_{\tau } \hat{\rho }_{L_{i-1}\bullet _i} \hat{O}^{i\dagger }_{\tau },\end{align}

\begin{matrix} {\hat{ρ}}_{L_{i - 1} •_{i}}^{'} = {\hat{ρ}}_{L_{i - 1} •_{i}} + a \sum_{τ} {\hat{O}}_{τ}^{i} {\hat{ρ}}_{L_{i - 1} •_{i}} {\hat{O}}_{τ}^{i †}, \end{matrix}

(14)

where a > 0 is an arbitrary constant, and

$\hat{O}^i_\tau$

{\hat{O}}_{τ}^{i}

is a site Hamiltonian operator (we use

$\hat{p}_i$

{\hat{p}}_{i}

⁠,

$\hat{p}_i^\dagger$

{\hat{p}}_{i}^{†}

⁠, and

$\hat{p}_i^\dagger \hat{q}_i$

{\hat{p}}_{i}^{†} {\hat{q}}_{i}

⁠). The constant a is used to control the strength of the perturbation during the sweep. We typically start with a value of 10⁻³ and decrease it by one order of magnitude when the energy from successive sweeps changes by less than 10⁻⁷E_h. When a reaches a minimum value (a < 10⁻⁷ 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.²² 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.²² Geometry optimization was performed using TURBOMOLE.²⁵ 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 space²⁶ indicate that ASD-DMRG with M = 128 recovered the CAS configuration interaction to within 10⁻⁷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¹⁶ total configurations, is intractable even for today's most powerful supercomputers.

FIG. 2.

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

View large Download slide

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

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 E_h).

		Δ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).²⁷ For monomers, we used geometries optimized using DFT under D_2h symmetry with the def2-TZVPP basis set²² 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⁷ and 2 × 10¹², 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 mE_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 mE_h. 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 E_h) 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 mE_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.³⁰

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

REFERENCES

1.

V.

Coropceanu

,

J.

Cornil

,

D. A.

da Silva Filho

,

Y.

Olivier

,

R.

Silbey

, and

J.-L.

Brédas

,

Chem. Rev.

107

,

926

(

2007

).

https://doi.org/10.1021/cr050140x

2.

M. B.

Smith

and

J.

Michl

,

Chem. Rev.

110

,

6891

(

2010

).

https://doi.org/10.1021/cr1002613

3.

S. M.

Parker

,

T.

Seideman

,

M. A.

Ratner

, and

T.

Shiozaki

,

J. Chem. Phys.

139

,

021108

(

2013

).

https://doi.org/10.1063/1.4813827

4.

S. M.

Parker

and

T.

Shiozaki

,

J. Chem. Theory Comput.

10

,

3738

(

2014

).

https://doi.org/10.1021/ct5004753

5.

J.

Olsen

,

B. O.

Roos

,

P.

Jørgensen

, and

H. J. Aa.

Jensen

,

J. Chem. Phys.

89

,

2185

(

1988

).

https://doi.org/10.1063/1.455063

6.

S. M.

Parker

,

T.

Seideman

,

M. A.

Ratner

, and

T.

Shiozaki

,

J. Phys. Chem. C

118

,

12700

(

2014

).

https://doi.org/10.1021/jp505082a

7.

S. R.

White

,

Phys. Rev. Lett.

69

,

2863

(

1992

).

https://doi.org/10.1103/PhysRevLett.69.2863

8.

S. R.

White

,

Phys. Rev. B

48

,

10345

(

1993

).

https://doi.org/10.1103/PhysRevB.48.10345

9.

S. R.

White

and

R. L.

Martin

,

J. Chem. Phys.

110

,

4127

(

1999

).

https://doi.org/10.1063/1.478295

10.

D.

Zgid

and

G. K.-L.

Chan

,

Annu. Rep. Comput. Chem.

5

,

149

(

2009

).

https://doi.org/10.1016/S1574-1400(09)00507-6

11.

Y.

Kurashige

and

T.

Yanai

,

J. Chem. Phys.

130

,

234114

(

2009

).

https://doi.org/10.1063/1.3152576

12.

K. H.

Marti

and

M.

Reiher

,

Z. Phys. Chem.

224

,

583

(

2010

).

https://doi.org/10.1524/zpch.2010.6125

13.

G. K.-L.

Chan

and

S.

Sharma

,

Annu. Rev. Phys. Chem.

62

,

465

(

2011

).

https://doi.org/10.1146/annurev-physchem-032210-103338

14.

S.

Wouters

and

D.

Van Neck

,

Eur. Phys. J. D

68

,

272

(

2014

).

https://doi.org/10.1140/epjd/e2014-50500-1

15.

Y.

Kurashige

,

G. K.-L.

Chan

, and

T.

Yanai

,

Nature Chem.

5

,

660

(

2013

).

https://doi.org/10.1038/nchem.1677

16.

S.

Sharma

,

K.

Sivalingam

,

F.

Neese

, and

G. K.-L.

Chan

,

Nature Chem.

6

,

927

(

2014

).

https://doi.org/10.1038/nchem.2041

17.

G.

Moritz

,

B. A.

Hess

, and

M.

Reiher

,

J. Chem. Phys.

122

,

024107

(

2005

).

https://doi.org/10.1063/1.1824891

18.

J.

Rissler

,

R. M.

Roack

, and

S. R.

White

,

Chem. Phys.

323

,

519

(

2006

).

https://doi.org/10.1016/j.chemphys.2005.10.018

19.

S. R.

White

,

Phys. Rev. B

72

,

180403

(

2005

).

https://doi.org/10.1103/PhysRevB.72.180403

20.

D.

Zgid

and

M.

Nooijen

,

J. Chem. Phys.

128

,

144116

(

2008

).

https://doi.org/10.1063/1.2883981

21.

J.

Pipek

and

P. G.

Mezey

,

J. Chem. Phys.

90

,

4916

(

1989

).

https://doi.org/10.1063/1.456588

22.

F.

Weigend

and

R.

Ahlrichs

,

Phys. Chem. Chem. Phys.

7

,

3297

(

2005

).

https://doi.org/10.1039/b508541a

23.

A. D.

Becke

,

J. Chem. Phys.

98

,

5648

(

1993

).

https://doi.org/10.1063/1.464913

24.

C.

Lee

,

W.

Yang

, and

R.

Parr

,

Phys. Rev. B

37

,

785

(

1988

).

https://doi.org/10.1103/PhysRevB.37.785

25.

F.

Furche

,

R.

Ahlrichs

,

C.

Hättig

,

W.

Klopper

,

M.

Sierka

, and

F.

Weigend

,

WIREs Comput. Mol. Sci.

4

,

91

(

2014

).

https://doi.org/10.1002/wcms.1162

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

).

https://doi.org/10.1016/j.tsf.2006.02.024

28.

D.

Zgid

and

M.

Nooijen

,

J. Chem. Phys.

128

,

014107

(

2008

).

https://doi.org/10.1063/1.2814150

29.

S.

Sharma

and

G. K.-L.

Chan

,

J. Chem. Phys.

136

,

144105

(

2012

).

https://doi.org/10.1063/1.3696963

30.

BAGEL, Brilliantly Advanced General Electronic-structure Library, See http://www.nubakery.org under the GNU General Public License.

© 2014 AIP Publishing LLC.

2014

AIP Publishing LLC