The *ab-initio* density matrix renormalization group (DMRG) is a tool that can be applied to a wide variety of interesting problems in quantum chemistry. Here, we examine the density matrix renormalization group from the vantage point of the quantum chemistry user. What kinds of problems is the DMRG well-suited to? What are the largest systems that can be treated at practical cost? What sort of accuracies can be obtained, and how do we reason about the computational difficulty in different molecules? By examining a diverse benchmark set of molecules: *π*-electron systems, benchmark main-group and transition metal dimers, and the Mn-oxo-salen and Fe-porphine organometallic compounds, we provide some answers to these questions, and show how the density matrix renormalization group is used in practice.

## I. INTRODUCTION

The density matrix renormalization group algorithm (DMRG), introduced by White,^{1,2} has by now been widely applied in quantum chemistry.^{3–15} Developing from its roots in condensed matter, its earliest application to chemical problems used the semi-empirical *π*-electron Pariser-Parr-Pople Hamiltonian.^{16–18} White and Martin introduced the first efficient formulation of the DMRG algorithm for *ab-initio* Hamiltonians (using some algorithmic contributions from Xiang^{19}) and the *ab-initio* density matrix renormalization group method was born. Since then, many groups have independently implemented and improved on the *ab-initio* DMRG algorithm. Some of these improvements include parallelization,^{8,20} non-Abelian symmetry and spin-adaptation,^{7,21–23} orbital ordering^{5,24–26} and optimization,^{9,27–29} more sophisticated initial guesses,^{5,24,25,30,31} better noise algorithms,^{5,32} extrapolation procedures,^{5,33,34} response theories,^{35,36} as well as the combination of the DMRG with various other quantum chemistry methods such as perturbation theory,^{37} canonical transformations,^{38} configuration interaction,^{39} and relativistic Hamiltonians.^{40}

In the ecosystem of quantum chemistry, the DMRG occupies a unique spot. On one hand, because the accuracy depends on a single, essentially continuously, tunable parameter—the number of renormalized states *M*—calculations can be converged and extrapolated in a simple way to the exact full configuration interaction limit. As a result, its earliest applications were as a proxy for full configuration interaction (FCI) in modest sized molecules where FCI would be too expensive. Examples include the early water^{34} and nitrogen^{41} calculations that used relatively small (but still inaccessible to FCI) basis sets, as well as more recent benchmarks, such as the exact treatment of the Be_{2} molecule^{42} at the complete basis set limit. On the other hand, the matrix product state (MPS) wavefunction^{43–46} underlying the DMRG is not based on an excitation expansion around a single-reference, and is thus well suited to non-mean-field, or strongly correlated, electronic structure, such as arising in transition-metal chemistry. Here, the DMRG has so far been applied in a complete active space setting,^{8,38,47–51} starting with the early work of Reiher and coworkers,^{47} to the latest calculations on systems as large as the bioinorganic Mn_{4}Ca core of photosystem II by Kurashige *et al.*,^{52} and the ubiquitous [4Fe-4S] biological iron-sulfur complexes by Sharma *et al.*;^{53} these have active spaces in excess of 50 orbitals. Finally, the internal structure of the MPS means that it is uniquely suited to pseudo-one-dimensional correlation, and in the *ab-initio* chemical context, the DMRG has been used since its inception to study ground- and excited-states of *π*-conjugated molecules,^{27,54,55} such as the polyacenes^{55} and graphene nanoribbons,^{49} as well as other one-dimensional systems, such as atomic chains and rings.^{54,56} The entanglement structure of the MPS has even generated a niche in interpretative quantum chemistry.^{52,57–60}

With the advent of publicly available quantum chemistry DMRG codes,^{10,20–22} the DMRG is transitioning into a method available not only to expert developers but also to the informed user. As with all complex methods, there is some degree of experience required to use it effectively. The purpose of this article is therefore to answer the following questions:

What sort of molecules can the DMRG be practically applied to?

What sort of accuracies can be obtained and at what cost? What are the typical sizes of systems (e.g., number of active orbitals) that can be treated with practical computational resources?

How do we reason about the accuracy of DMRG calculations for different molecules?

How is a DMRG calculation best specified (e.g., in terms of starting orbitals and their order)?

We answer these questions both from theoretical reasoning and by applying the method to a “representative” set of molecules for DMRG applications: *π*-conjugated systems (polyenes and arenes), benchmark dimers (C_{2} and Cr_{2}) in several basis sets (with more than 100 orbitals), and organometallic complexes (Fe-porphine and a Mn-oxo-salen model of Jacobsen’s catalyst) in large active spaces (with more than 40 active orbitals). Importantly, the calculations we describe are all run in a completely black-box fashion using the default settings of our Block code, and thus should be easily replicated. The Block code is currently available in the Molpro,^{61}^{,}Orca,^{62} and Q-Chem^{63} distributions, and can also be freely downloaded from our weblink.^{64}

In Secs. II A–II E, we give a self-contained overview of the DMRG algorithm and MPS. We describe both the theory and how to reason about it with respect to orbital choices, ordering, accuracy, and cost, from the viewpoint of the user. We further define some of the choices and settings used in our particular implementation of the DMRG algorithm in Block, with respect to the sweep and noise schedule. The various issues are then explored in the context of the set of “representative” molecules in Sec. III. The summary of our findings and recommendations is given in Sec. IV.

## II. THEORETICAL BACKGROUND

### A. The DMRG wavefunction and optimization algorithm

As with other wavefunction methods in chemistry, the DMRG is based on an approximate wavefunction ansatz. This ansatz is known as the MPS. A MPS is a non-linear wavefunction ansatz, built as a product of variational objects for each orbital in the basis. Technically, the DMRG refers to the combination of a specific, efficient, and self-consistent optimization algorithm (the sweep algorithm) with the variational MPS ansatz. In fact, historically, the DMRG is discussed primarily in terms of the optimization algorithm,^{1,2} which derives from the original numerical renormalization group,^{65} rather than in terms of its wavefunction structure. However, the two pictures are closely intertwined because in practice the existence of an efficient optimization scheme for the highly non-linear ansatz is a key to the success of the DMRG in practice. The close relationship between the wavefunction and optimization method in DMRG recalls the close relationship between the Hartree-Fock Slater determinant approximation and the self-consistent field algorithm used to optimize it. Indeed, in quantum chemistry, we often interpret “SCF” to mean Slater determinant in an eigenvalue (e.g., ground-state) context, and similarly “DMRG” is often synonymous with MPS for eigenvalue problems. Note that the parallel terminology is not the only similarity between DMRG and Hartree-Fock theory; the product structure of the MPS is the source of further useful analogies, discussed in detail in Refs. 66–69.

In practice, a DMRG calculation will often use two types of matrix product states: a *one-site* and a *two-site* MPS. (Site is here synonymous with orbital, in an occupation number representation.) These refer to different but related ansatz. The one-site MPS is easier to reason about and converge to high precision, but optimization with this ansatz alone often ends up trapped in local minima^{32} due to the inability to change “quantum numbers” (a type of internal wavefunction symmetry) during the DMRG sweep, as described in Sec. II C. A standard strategy is therefore to start the DMRG calculation using the more flexible two-site MPS and then to switch near convergence to the one-site MPS.

To define a MPS, we are first required to order the orbitals *ϕ*_{1}…*ϕ _{k}*. Then, the one-site MPS is defined as

Here, *n _{i}* denotes the occupancy of orbital

*ϕ*. It takes four values corresponding to the states | − 〉, |↑〉, |↓〉, and |↑↓〉. For a given index

_{i}*n*,

_{i}**A**

^{ni}is then an

*M*×

*M*matrix, and each

**A**is correspondingly a 3-index tensor.

**A**

^{ni}contains information on the wavefunction coefficients for that particular choice of occupancy of orbital

*ϕ*. The first and last

_{i}**A**’s are slightly special, as the dimensions of

**A**

^{n1}and

**A**

^{nk}are 1 ×

*M*and

*M*× 1, respectively; this ensures that for a given occupancy string

*n*

_{1}…

*n*, performing the vector, matrix, matrix …vector product defined by

_{k}**A**

^{n1}

**A**

^{n2}…

**A**

^{nk}yields a scalar wavefunction amplitude, corresponding to the coefficient of the determinant |

*n*

_{1}…

*n*〉.

_{k}We see that the DMRG is a kind of product wavefunction, but unlike a Slater determinant, there is a variational object associated with each *orbital in the basis*, rather than with each electron. The variational freedom and thus accuracy of a DMRG wavefunction is controlled by the single dimensional parameter *M*, referred to as the “number of renormalized states,” or MPS bond dimension. In total, there are $O ( M 2 k ) $ variational parameters in the wavefunction. When *M*^{2} ∼ *N _{d}* (where

*N*is the total number of possible determinants), the DMRG wavefunction becomes exact, but typically the convergence with

_{d}*M*is rapid and a much smaller

*M*(in the range of 1000–100 00) is used for chemical accuracy in practice.

The two-site DMRG wavefunction is defined similarly, but the **A** tensors of a special pair of adjacent (“two”) sites are fused together to become a four-index tensor **T**

For each pair of occupancies *n _{i}* and

*n*

_{i+1},

**T**

^{nini+1}is an

*M*×

*M*matrix. The

**T**tensor in the two-site MPS gives it slightly more variational freedom than the one-site MPS. This allows for an optimization of the quantum numbers as described below, which helps to avoid local minima. During the DMRG sweep the boundary for the special pair of sites is iterated from the first (last) pair of orbitals to the last (first) during the DMRG forwards (backwards) sweep optimization; each boundary corresponds to a single step of the sweep. Thus in the DMRG sweep there is not a unique two-site MPS, but rather

*k*− 1 of them, depending on where we choose to place the two sites

*i*+ 1. When we refer to the DMRG energy from a two-site wavefunction, we therefore mean the energies obtained at the lowest point in the sweep. This is often near to when the two sites are at the middle of the lattice.

### B. Sweeps and operators

In the DMRG sweep algorithm, the energy is variationally optimized with respect to a *single* **A** (or **T** in the case of the two-site MPS) tensor at a time. This is analogous to Hartree’s original SCF procedure in Hartree-Fock theory where a single orbital is updated at a time, and is also conceptually related to “alternating” methods in multi-linear algebra. The sweep traverses tensors from left to right (forwards) and from right to left (backwards) and multiple sweeps are carried out to converge the energy. To obtain the updated **A**^{ni} matrices at each site, we solve an effective Schrödinger equation with an effective Hamiltonian (analogous to a Fock equation for each site) for the $O ( M 2 ) $ coefficients in **A**^{ni}. Because *M*^{2} can be quite large, and we only need the lowest eigenvalue and eigenvector (if we are interested in a single state), we use the Davidson algorithm.^{70} The wavefunction solution is the dominant computational cost in a DMRG calculation. It has a formal cost scaling of $O ( M 3 k 3 ) $ per sweep, where *M*^{3} derives from matrix multiplication, although in practice $O ( M 2 ) $ scaling is often observed due to the block-sparse structure of the *A*^{ni} matrices deriving from quantum numbers, as is used in this work.

The effective Hamiltonian at each site is defined by tracing the Hamiltonian in the energy expectation value 〈Ψ|*H*|Ψ〉 with tensors from all other sites in the bra and ket. The efficient construction of the effective Hamiltonian requires decomposing it into contributions from all orbitals to the left of the current optimized tensor (the left block), and orbitals to the right of the current optimized tensor (the right block). These intermediates (renormalized operators) are stored in memory and on disk, and form the dominant memory cost of the DMRG calculation: $O ( M 2 k 2 ) $ memory and $O ( M 2 k 3 ) $ disk.

Because the sweeps are variational in nature, the DMRG energy at each sweep decreases monotonically. (Strictly speaking, monotonicity is guaranteed only with the one-site algorithm but is usually observed with the two-site algorithm also except when very close to final convergence.) It is more efficient to carry out the initial DMRG sweeps with small values of *M* to approximately converge the DMRG wavefunction, and then to carry out later sweeps with larger *M*. The choice of successive increasing *M*’s leads to a *sweep schedule*. The default sweep schedule in Block, as used in all our calculations here, is listed in Table I. In the calculation, the desired final *M* is input, and the sweep schedule is followed up to the desired *M*. The first sweep (the warmup sweep) is a little bit special as we need to provide an initial guess for the **A** tensors in the MPS. In Block, we build the **A** tensors for the warmup sweep corresponding to a small number of low-energy determinants.^{5}

M . | No. of iterations . | Noise . |
---|---|---|

500 | 8 | 1 × 10^{−4} |

1000 | 8 | 5 × 10^{−5} |

2000 | 4 | 5 × 10^{−5} |

3000 | 4 | 5 × 10^{−5} |

… | 4 | 5 × 10^{−5} |

Max. M | 4 | 5 × 10^{−5} |

Max. M | 2 | 0 |

M . | No. of iterations . | Noise . |
---|---|---|

500 | 8 | 1 × 10^{−4} |

1000 | 8 | 5 × 10^{−5} |

2000 | 4 | 5 × 10^{−5} |

3000 | 4 | 5 × 10^{−5} |

… | 4 | 5 × 10^{−5} |

Max. M | 4 | 5 × 10^{−5} |

Max. M | 2 | 0 |

### C. Quantum numbers, symmetries, and noise

Global symmetries, including Abelian symmetries (such as particle number, point-group symmetry, and *S _{z}* symmetry) as well as non-Abelian symmetries (such as the SU(2) or

*S*

^{2}symmetry, or non-Abelian point-group symmetry) can be efficiently used with the DMRG algorithm and our Block implementation. We can ensure that the underlying MPS transforms according to an irrep of a global symmetry, by assigning irreducible representations to the left and right matrix indices of

**A**

^{ni}. In this context, the irreps associated with the left and right matrix indices are referred to as “quantum numbers.” The global irrep then implies a block-sparse structure in the

**A**

^{ni}matrices. The global molecular symmetry thus imply symmetries of local transformations of the matrices.

Note that symmetry does not itself dictate the *sizes* of the non-zero blocks of **A**^{ni}. In principle, these sizes should be optimized during the optimization of the wavefunction, as an additional discrete optimization problem. A poor distribution of block sizes in the tensors is the main cause of the local minimum problem in DMRG, where a calculation appears to converge to a state of too high energy. This phenomenon is sometimes also referred to as having incorrect, or “losing” quantum numbers (i.e., not having the appropriate irrep labels for the states).^{20} One way to detect an incorrect convergence to a local minimum is to increase *M*. A calculation with incorrect quantum numbers at small *M* will exhibit a very sudden lowering of the energy at larger *M* when the correct quantum number sector is finally recovered.

So far in the literature, the discrete optimization of block sizes has not been addressed in a very sophisticated way in DMRG algorithms. Since the original formulation it was recognized that the two-site MPS allowed for a limited local variation of block sizes for the particular choice of the “two sites.”^{2,5} When coupled with perturbative noise which introduces random quantum numbers into the wavefunction, a DMRG sweep with the two-site wavefunction provides some ability to dynamically search different block distributions.^{32} Using the perturbative noise is particularly important in the initial sweeps in the DMRG, as different block distributions can be easily ignored due to the small total number of states *M*. However, the noise should be gradually reduced in later sweeps in the schedule to avoid affecting the final converged answer. In conjunction with the default sweep schedule in Block, we use an accompanying noise schedule as shown in Table I. For the final sweeps at the desired maximum *M*, the noise is set to zero.

### D. Truncation error and extrapolation

Because the two-site MPS has additional variational freedom over the one-site MPS, there is, in general, an error from projecting from a two-site MPS to a one-site MPS, which only vanishes in the limit of large *M*. This “truncation error” (equivalent to the discarded weight in the density matrix in the traditional DMRG language) provides an estimate of the accuracy of the DMRG wavefunction. It was empirically observed^{5,33} that the truncation error is almost linearly proportional to the error in the DMRG energy. In an exact (i.e., for sufficiently large *M*) calculation, the truncation error becomes zero, and thus extrapolating to zero error gives an estimate of the “exact” energy for infinite *M*.

As a practical note, the “noise” in the initial sweeps introduces random error into the DMRG wavefunction and can cause problems with extrapolation when using the data at small *M*. One can avoid this by carrying out a few sweeps at each *M* without noise before proceeding to higher *M*.^{10} Alternatively, after the calculation is converged, we can carry out a “reverse schedule” where *M* starts at its maximum value and is then *lowered* to its smallest value, to obtain good data at small *M* (which is often biased by the warmup sweep). We call this the *reverse schedule*. Using energies from the reverse schedule can sometimes lead to more accurate linear extrapolations of the DMRG energy. Comparison of the extrapolations using the standard schedule and reverse schedule data is shown in Fig. 1.

### E. Orbital choice and ordering in DMRG

Choosing an orbital space for a DMRG calculation, much like the selection of an active space in other wavefunction methods, requires some trial and error. However, unlike many other common wavefunction methods, a MPS is not invariant to orbital rotations within the active space, except at large *M*.^{27,28} Orbital re-ordering can be viewed as a large rotation. Although the DMRG calculation becomes less sensitive to these choices as *M* increases, at fixed *M*, a good choice of orbitals and ordering greatly improves the accuracy of a DMRG calculation in practice.

In principle, the best ordering and choice of orbitals minimizes the maximum entanglement at any cut of the DMRG orbital lattice (the lattice being the one-dimensional ordering of orbitals). Indeed, for a MPS wavefunction with *M* renormalized states, there is an upper bound of log_{2} *M* on the amount of entanglement that can be captured by such a state. A related criterion is to minimize entanglement (defined in some way) between distant orbitals in the lattice.

Entanglement is distinct from quantum chemical correlation. One way to reason about ground-state entanglement is to view it as generated by the one- and two-electron parts of the Hamiltonian, as we project from a simple product state of orbitals to the ground-state with the operator exp(−*βH*), *β* → ∞. If we start with very localized orbitals such as orthogonalized atomic orbitals, the one-electron part of the Hamiltonian will cause the orbitals to mix due to electron delocalization (to reduce the “spread” in the energy domain) and this can generate a large amount of entanglement. The “one-electron” entanglement effects can be removed by transforming to molecular orbitals. On the other hand, the two-electron part of the Hamiltonian generates entanglement also, which can typically be minimized by working with a more local (in real-space) representation, as familiar from local correlation methods. A reasonable compromise is to use split-localized orbitals, which correspond to orbitals where the occupied and virtual orbitals are separately localized.

Once the orbitals have been determined, we still need to choose an order. (One could, in principle, optimize both the shapes and the order of the orbitals directly via energy minimization: in complete active space terminology, this would correspond to “active-active” rotations. We do not consider these in this work.) The ordering should place orbitals which are most entangled (such as orbitals which are close in space, or pairs of bonding and anti-bonding orbitals) close together. In quasi-one-dimensional molecules (such as chain, strip, or ring-like molecules), there is an obvious ordering that follows the connectivity of the molecule. However, finding the optimal ordering in general is NP-hard.

Nonetheless, reasonable heuristic algorithms are available. These start with by defining a matrix, or *entanglement metric*. A rigorous measure of entanglement between two orbitals is the mutual information (pair entanglement entropy) first investigated by Rissler and White.^{26} However, this requires knowledge of the final wavefunction. A much simpler proxy, which has some of the correct properties, is the exchange integral between the two orbitals $ K i j =\u222b\u222bd r 1 d r 2 r 12 \u2212 1 \varphi i \u2217 ( r 1 ) \varphi j ( r 2 ) \varphi j \u2217 ( r 1 ) \varphi i ( r 2 ) $, which measures their proximity and spatial overlap. (The one-electron Hamiltonian was initially used in Ref. 5 and provides an alternative simple metric.) The goal is then to optimize the ordering such that $ \u2211 i j | K i j |$ (or some related function) is minimized.

In our Block code, we have two default orderings: a Fiedler vector order and a genetic algorithm optimized order. The Fiedler vector is a graph theoretic technique which provides good approximations to the graph min-cut problem, which is what we are interested in.^{71–73} Graph techniques had earlier been examined in the DMRG ordering problem,^{5} but a detailed study of the Fiedler vector was first presented by Barcza *et al.* in Ref. 25. The derived ordering can be obtained as follows. First, we assume that the positions of the *k* orbitals on the DMRG orbital lattice are continuous 1D variables *x*_{1}…*x _{k}*. Then, we measure the distance between orbitals on the lattice by

We now minimize *F*(*x*) with respect to the vector of positions {*x*}. To prevent rigid translation of all “coordinates,” we fix $ \u2211 i x i =0$, and to prevent a trivial solution where *x* = 0, we impose normalization $ \u2211 i x i 2 =1$. We therefore solve min *F*({*x*}) subject to $ \u2211 i x i =0$; $ \u2211 i x i 2 =1$, a simple linear algebra problem, equivalent to finding the second lowest eigenvector of the graph Laplacian **L** = **D** − **K** where **D** is a diagonal matrix with entries $ D i i = \u2211 j | K i j |$ (the first eigenvector is simply the constraint condition).^{74,75} Sorting the values of the vector coordinates {*x*} then gives the Fiedler ordering. The advantage of this method is that it is trivial to obtain from a small matrix diagonalization, and since our optimization metric |*K _{ij}*| is in any case approximate, it is not necessary to obtain the true global optimum. For this reason, the Fiedler ordering is the default ordering in the Block code.

We can look for additional improvement in the metric by using a genetic algorithm (GA). Genetic algorithms are global optimization algorithms based on an analogy to natural selection^{76} and were first employed in the DMRG orbital ordering problem by Moritz *et al.*^{77} In our work, we use as our cost function

where *D _{ij}*(

*x*) is the distance matrix that depends on the ordering

*x*. To implement the GA, we need crossover and/or mutation operators upon selection. The crossover operator generates a new ordering from two different input ordering strings. The implementation in Block uses a

*partially matched crossover*(PMX), which works well for permutation strings.

^{78,79}The mutation operator modifies a string randomly, and in our case, we pick two parts of string randomly and swap them to create a new string. We have further chosen the probabilities of selection (by trial and error) to avoid local minima for the orbital ordering problem. We start the GA using the Fiedler vector as an initial guess, which greatly reduces the optimization time in the GA. For all cases we tested, the solution with and without the Fiedler vector guess results in the same GA minimum.

### F. Accuracy and computational cost

As discussed above, once the orbitals and ordering are defined, the accuracy of the DMRG calculation is controlled by a single parameter *M*, the number of renormalized states. The memory and CPU costs are summarized as follows:

Memory: The scaling with the number of states

*M*and orbitals*k*is roughly*O*(*M*^{2}*k*^{2}), while disk usage is*O*(*M*^{2}*k*^{3}). DMRG calculations are very memory intensive. The largest calculation in this work is the butadiene calculation (22*e*, 82 orbs,*C*_{1}symmetry), which used 850 Gb of memory and 17 Tb of disk. However, note that the parallel DMRG algorithm used in the Block code^{20}distributes both memory and disk storage across different nodes, and thus even modest sized computational clusters provide quite substantial resources.CPU: DMRG CPU times scale as

*O*(*M*^{3}*k*^{3}) per sweep. For small*M*(say, less than 1000), and for the initial DMRG warmup sweep, our implementation has a large*M*^{2}overhead associated with large sets of quantum numbers which can dominate the computational cost. Using the parallelization strategy in Ref. 20, we observe reasonable parallelization up to a number of cores equal to the number of orbitals. For the largest calculation in this work (butadiene, 22*e*, 82 orbs,*C*_{1}symmetry), a single sweep at*M*= 3000 (which gave an energy below CCSDT) took 25 h on 42 cores.

In almost all cases in chemistry, the DMRG energy converges almost exponentially with *M* (empirically, a convergence of exp(*κ*(ln*M*)^{α}), where *α* ≈ 2 is observed^{5,80}); however, the exponent *κ* of the exponential is molecule dependent. To understand *κ* and the rate of convergence in different molecules, we can make a few general statements. As discussed above, *M* determines the maximum amount of entanglement by the MPS wavefunction. In the simple case of (non-critical) 1D system where every orbital is equivalent (e.g., a chain of hydrogen atoms in a minimal orthogonalized atomic basis) the ground-state entanglement is controlled only by the gap of the system and is otherwise independent of length. In this case, *M* *can be held fixed and will give a fixed accuracy per electron* regardless of system size, and the only scaling of the method derives from the significant Hamiltonian matrix elements. This is why the DMRG and MPS are informally said to give a polynomial complexity solution for correlation in one-dimensional systems. Moving to two-dimensional systems in a similar atomic basis (i.e., where every orbital is equivalent), the amount of entanglement grows *linearly* with the width of the system (for systems that satisfy the area law). This requires *M* to depend *exponentially* on width, which is why the DMRG is said to have an exponential scaling with the width of a 2D system.

However, in many molecular calculations, we may be more interested in the computational scaling with respect to increasing the basis size (for fixed energy accuracy and molecular size). In this case, as one is increasing the basis, the orbitals entering are not all equivalent; for example, one might have in the basis valence orbitals and correlating Rydberg like orbitals. These typically have very different occupancies and contributions to the wavefunction and energy. In a large basis, the molecular wavefunction, even in strongly correlated systems, tends to be somewhat “concentrated” around the occupied and valence orbitals. Linear combinations of a few (say *N*_{large}) large weighted determinants, as is the case for many multireference problems in molecules, can only generate log*N*_{large} entanglement and are thus easy to capture with a MPS. It can be more effort to capture the manifold of double excitations that, in small molecules, provide the bulk of the dynamical correlation contributions. Assuming that the reordering of orbitals results in the occupied and virtual orbitals being distributed randomly in the DMRG lattice, this generates $O ( log n occ n virt ) $ entanglement, and thus *M* ∝ *n*_{occ}*n*_{virt}. Thus, we expect a worst case scaling of $O ( n virt 6 ) $, for fixed accuracy in the dynamical correlation as the basis size increases, for a fixed molecular size. (As the molecular size increases, the locality results in a lower scaling, as reasoned about in the paragraphs above.)

## III. RESULTS

### A. Arenes

We begin with the arenes, shown in Fig. 2. These conjugated molecules can be viewed as finite strips of graphene of various “widths.” We used octatetraene (8 carbon atoms), tetracene (18 carbon atoms), and three larger analogous poly-aromatic hydrocarbons with 28, 38, and 48 carbon atoms, respectively, corresponding to widths **1**–**5**. The range of widths allows us to study the dependence of the DMRG energy convergence on molecular dimensionality. To provide an idealized test system for orbital shapes and orbital ordering, we limit this particular study to the *π*-electron system in a minimal (STO-3G) basis.^{81} We used optimized restricted Hartree-Fock (RHF/STO-3G) geometries. Within the *π*-active spaces, corresponding to (8*e*, 8*o*), (18*e*, 18*o*), (28*e*, 28*o*), (38*e*, 38*o*), and (48*e*, 48*o*) for **1**–**5**, respectively, we considered several orbital choices: canonical Hartree-Fock orbitals, split-localized orbitals (that resemble two-centre *π*-bonding and anti-bonding orbitals), and fully-localized orbitals (orthogonalized *p* atomic orbitals). We localized orbitals using the Pipek-Mezey method.^{82}

For each orbital choice, the subsequent ordering was determined using both the Fiedler vector and GA optimization schemes described in Sec. II E. These orderings are shown for localized orbitals in Fig. 3. It is known that a good DMRG ordering in 2D sheets is a simple “snake-like” pseudo-1D ordering traversing each row.^{83} This is reproduced in the smaller arenes both by the Fiedler and GA orderings. For arene **5**, which has more of a square aspect ratio, the Fiedler and GA orderings produce orderings that roughly go diagonally across the arene, and the resulting orderings are slightly different.

Fig. 4 shows the energy convergence for arene **3** for different orbital choices and orderings. This convergence is representative of the other arenes. The differences between the Fiedler and GA ordering results are very small; the effect of orbital choice is far more significant. Due to the valence only nature of the active space, we observe fastest convergence with fully-localized orbitals (atomic-like *p* orbitals), even though the energy at small *M* is relatively poor. In contrast, the canonical Hartree-Fock orbitals give good energies at very small *M* (since the Hartree-Fock limit is captured by *M* = 1) but only slowly converge the energy as *M* is increased. Split-localized orbitals are a good compromise, capturing the Hartree-Fock limit at *M* = 1 and giving better energies than the localized orbitals at very small *M*, while giving rapid convergence to the exact result that is only slightly slower than the fully localized basis.

Note in Fig. 4 the near-exponential convergence of the energy error with *M*, typical of the DMRG. As discussed in Sec. II F, a more careful analysis predicts that the energy error follows $ln \delta E \u223c\kappa ln M 2 $.^{5} This behaviour is seen from the plot in Fig. 5. Finally, the linear relationship between the discarded weight and the energy is shown in Fig. 6. This demonstrates the simple extrapolation of the energy to zero discarded weight.

In Fig. 7, we study the DMRG energy convergence as a function of the arene “width.” As expected, the *M* required for a given accuracy grows roughly exponentially with width, reflecting the exponentially increasing set of fluctuations across any “boundary” that cuts across the system. However even for the largest arene, we see that it is possible to converge the active space energy per carbon atom in 48 orbitals nearly exactly (to 0.01 m*E _{h}*) with an accessible

*M*< 100 00. (As a technical note, it should be noted here that because we use

*S*

^{2}symmetry in our DMRG code, our

*M*values refer to the number of multiplets retained. The above would thus correspond to

*M*= 200 00 − 300 00 states in a code with only

*S*symmetry.)

_{z}### B. Dimers

In our next test case, we consider two diatomics, C_{2} and Cr_{2}. Diatomics provide a contrast to the *π*-electron system of the arenes as there is little meaningful spatial locality in such a small molecule. Further here, we carry out calculations with all electrons and larger basis sets. This provides a wide range of energies and weights in the wavefunction, unlike in the arenes, and allow us to explore the use of the DMRG for dynamical correlation.

We start with C_{2}. We use the equilibrium bond distance of 1.24253 Å (the same as in the initiator full configuration interaction quantum Monte Carlo (*i*-FCIQMC) calculations in Ref. 84) and perform calculations in Dunning’s cc-pVDZ, cc-pVTZ basis, and cc-pVQZ basis sets. We carry out both all electron (AE) calculations (corresponding to (12*e*, 28*o*), (12*e*, 60*o*), and (12*e*, 110*o*) spaces for cc-pVDZ, cc-pVTZ, and cc-pVQZ, respectively). We first examine the effect of different orderings. There is no natural ordering in a diatomic molecule, unlike for localized orbitals in the arenes. However, we can still reorder the orbitals using the Fiedler and genetic algorithms as for the arenes. The corresponding energies are shown in Table II. Also shown in Table II are energies obtained using the reverse sweep schedule where *M* is decreased from its maximum value of 5000 to 500. These are much better at smaller *M* than the energies obtained in the standard sweep, evidence that the warmup sweep is relatively poor in the standard schedule and also that additional important quantum numbers are detected during the sweeps at larger *M* that were not picked at smaller *M*.

. | Energy . | Genetic . | Fiedler . | Reverse genetic . |
---|---|---|---|---|

M . | E (E_{h})
. | E (E_{h})
. | E (E_{h})
. | E (E_{h})
. |

500 | −0.846 422 | −0.839 705 | −0.838 761 | −0.853 141 |

1000 | −0.851 769 | −0.847 877 | −0.849 565 | −0.855 803 |

2000 | −0.854 534 | −0.853 348 | −0.853 724 | −0.856 866 |

3000 | −0.855 732 | −0.855 104 | −0.855 446 | −8.857 135 |

4000 | −0.856 425 | −0.855 890 | −0.856 238 | −0.857 231 |

5000 | −0.856 762 | −0.856 386 | −0.856 645 | −0.857 256 |

Extrapolated | −0.857 28 |

. | Energy . | Genetic . | Fiedler . | Reverse genetic . |
---|---|---|---|---|

M . | E (E_{h})
. | E (E_{h})
. | E (E_{h})
. | E (E_{h})
. |

500 | −0.846 422 | −0.839 705 | −0.838 761 | −0.853 141 |

1000 | −0.851 769 | −0.847 877 | −0.849 565 | −0.855 803 |

2000 | −0.854 534 | −0.853 348 | −0.853 724 | −0.856 866 |

3000 | −0.855 732 | −0.855 104 | −0.855 446 | −8.857 135 |

4000 | −0.856 425 | −0.855 890 | −0.856 238 | −0.857 231 |

5000 | −0.856 762 | −0.856 386 | −0.856 645 | −0.857 256 |

Extrapolated | −0.857 28 |

Note also that the energy at *M* = 5000 for the reverse schedule (GA ordering) is about 0.9 m*E _{h}* lower than the

*M*= 5000 energy (GA ordering) in the default schedule. This difference is entirely due to the noise in the default sweeps, which is turned off for the reverse schedule. (Recall that noise is used to ensure more robust convergence in the early sweeps.) We see that in this molecule, the default setting for the noise in Table I should be decreased at larger

*M*to achieve faster convergence. This is unsurprising as the sweep parameters have been chosen for general utility rather than for specific molecules and illustrates the gains in convergence that could in principle be made by optimizing the sweep settings on a per molecule basis. However, the detrimental effect of the non-optimal noise remains quite small and is unimportant on the chemical scale of 1 m

*E*.

_{h}For the GA ordering, we further give converged energies (accurate to better than 0.05 m*E _{h}*) for the 1

*s*frozen core (FC) calculations in the same bases. This allows direct comparison (shown in Table III) to

*i*-FCIQMC calculations in the literature.

^{84}The

*i*-FCIQMC energies reported in Ref. 84 for the TZ and QZ bases are in rough agreement with the DMRG energies. However, the variational DMRG energies are clearly lower and have significantly smaller error bounds. The difference between the DMRG and

*i*-FCIQMC energy corresponds to the remaining initiator error in the

*i*-FCIQMC calculation, as we have also found in other comparisons between the two techniques.

^{42}

M . | DZ . | TZ . | QZ . |
---|---|---|---|

500 | −0.731 449 | −0.807 296 | −0.853 141 |

1000 | −0.731 856 | −0.808 662 | −0.855 803 |

2000 | −0.731 945 | −0.809 123 | −0.856 866 |

4000 | −0.731 958 | −0.809 260 | −0.857 231 |

6000 | … | −0.809 285 | … |

(FC) DMRG | −0.728 556 | −0.785 054 | −0.802 671 |

(FC) i-FCIQMC | −0.7287(8) | −0.7849(3) | −0.8025(1) |

M . | DZ . | TZ . | QZ . |
---|---|---|---|

500 | −0.731 449 | −0.807 296 | −0.853 141 |

1000 | −0.731 856 | −0.808 662 | −0.855 803 |

2000 | −0.731 945 | −0.809 123 | −0.856 866 |

4000 | −0.731 958 | −0.809 260 | −0.857 231 |

6000 | … | −0.809 285 | … |

(FC) DMRG | −0.728 556 | −0.785 054 | −0.802 671 |

(FC) i-FCIQMC | −0.7287(8) | −0.7849(3) | −0.8025(1) |

In C_{2}, we are using basis sets with many virtual orbitals, thus much of the correlation we describe is dynamic. This, together with the lack of spatial locality in the system, means (as described in Sec. II F) that the *M* required for a given accuracy in a molecule should scale roughly as *n _{occ}n_{virt}* to capture the dominant doubles excitations, i.e., the required

*M*should scale linearly with the number of electrons and basis size. The

*M*required to achieve an accuracy of 0.01 m

*E*is plotted against

_{h}*n*in Fig. 8. We observe a good linear-fit to

_{occ}n_{virt}*n*, confirming that the entanglement is dominated by the double excitations, and that correlation in this molecule is mainly dynamic in character.

_{occ}n_{virt}We next consider the Cr_{2} dimer. Cr_{2}, with its formal hextuple bond, is a long-standing favourite of quantum chemistry, due to the difficulty in adequately describing the shelf-like nature of its potential energy curve.^{85} This unusual curve arises from a combination of non-dynamic (strong) correlation and dynamic correlation in large basis sets. Cr_{2} has been studied previously with DMRG methods, which (in combination with perturbation theory) give a very accurate description of the Cr_{2} potential energy curve.^{37}

Here, for benchmark purposes, we use smaller basis sets where it is possible to converge the DMRG energy “exactly,” i.e., to beyond chemical accuracy. We consider Cr_{2} at the Cr-Cr distance of 1.5 Å with the Ahlrichs’ SV basis,^{86} both correlating all electrons (a (48*e*, 42*o*) space) as well as a (24*e*, 30*o*) active space subset, used in previous DMRG benchmarks.^{21,37} We used GA ordering for all calculations.

The DMRG energies and weights are shown in Table IV with CCSD(T) and CCSDTQ calculations for comparison. We see a large difference between the CCSDTQ and DMRG energy indicating the importance of high-order correlations in this complex system. The linear extrapolation of the energies in Table IV is plotted in Figure 9 from which obtain the estimated exact energies. The error bars from the fitting are ±0.007 and ±0.19 m*E _{h}* for the (24

*e*, 30

*o*) and (48

*e*, 42

*o*) active spaces, respectively, probably an underestimate. Alternatively, a conservative rule of thumb used in some DMRG extrapolations

^{87}is to assign the error bar as 1/5 the extrapolation energy, which comes to ±0.034 and ±0.32 m

*E*for the two active spaces, respectively, probably overestimates. Either way, the energies are very accurate. Further, with the maximum

_{h}*M*= 8000, even the

*unextrapolated*energies are within 0.2 m

*E*of the estimated exact result in the (24

_{h}*e*, 30

*o*) active space, and within 1.6 m

*E*of the estimated exact result in the (48

_{h}*e*, 42

*o*) active space.

. | (24e, 30o)
. | (48e, 42o)
. | ||||||
---|---|---|---|---|---|---|---|---|

. | Default . | Reverse . | Default . | Reverse . | ||||

M . | DW . | E (E_{h})
. | DW . | E (E_{h})
. | DW . | E (E_{h})
. | DW . | E (E_{h})
. |

500 | 1.02 × 10^{−4} | −0.414 935 | 9.47 × 10^{−5} | −0.416 106 | 2.73 × 10^{−4} | −0.177 377 | 4.78 × 10^{−4} | −0.421 928 |

1000 | 4.59 × 10^{−5} | −0.418 219 | 4.36 × 10^{−5} | −0.418 693 | 2.18 × 10^{−4} | −0.426 115 | 2.89 × 10^{−4} | −0.432 057 |

2000 | 2.34 × 10^{−5} | −0.419 707 | 1.87 × 10^{−5} | −0.419 986 | 1.34 × 10^{−4} | −0.436 820 | 1.51 × 10^{−4} | −0.438 127 |

3000 | 1.38 × 10^{−5} | −0.420 246 | 1.10 × 10^{−5} | −0.420 386 | 9.24 × 10^{−5} | −0.439 855 | 9.59 × 10^{−5} | −0.440 324 |

4000 | 1.10 × 10^{−5} | −0.420 480 | 7.78 × 10^{−6} | −0.420 570 | 7.14 × 10^{−5} | −0.441 211 | 6.92 × 10^{−5} | −0.441 458 |

5000 | 9.26 × 10^{−6} | −0.420 609 | 5.65 × 10^{−6} | −0.420 672 | 5.55 × 10^{−5} | −0.441 927 | 5.41 × 10^{−5} | −0.442 146 |

6000 | 6.69 × 10^{−6} | −0.420 686 | 4.44 × 10^{−6} | −0.420 735 | 4.50 × 10^{−5} | −0.442 439 | 4.37 × 10^{−5} | −0.442 607 |

7000 | 5.58 × 10^{−6} | −0.420 745 | 3.44 × 10^{−6} | −0.420 776 | 3.61 × 10^{−5} | −0.442 792 | 3.66 × 10^{−5} | −0.442 933 |

8000 | 4.55 × 10^{−6} | −0.420 774 | 2.69 × 10^{−6} | −0.420 780 | 3.15 × 10^{−5} | −0.443 334 | 3.08 × 10^{−5} | −0.443 173 |

Extrapolated | 0.0 | −0.420 948 | 0.0 | −0.444 784 | ||||

CCSD(T) | −0.398 638 | −0.422 229 | ||||||

CCSDTQ | −0.406 696 | −0.430 244 |

. | (24e, 30o)
. | (48e, 42o)
. | ||||||
---|---|---|---|---|---|---|---|---|

. | Default . | Reverse . | Default . | Reverse . | ||||

M . | DW . | E (E_{h})
. | DW . | E (E_{h})
. | DW . | E (E_{h})
. | DW . | E (E_{h})
. |

500 | 1.02 × 10^{−4} | −0.414 935 | 9.47 × 10^{−5} | −0.416 106 | 2.73 × 10^{−4} | −0.177 377 | 4.78 × 10^{−4} | −0.421 928 |

1000 | 4.59 × 10^{−5} | −0.418 219 | 4.36 × 10^{−5} | −0.418 693 | 2.18 × 10^{−4} | −0.426 115 | 2.89 × 10^{−4} | −0.432 057 |

2000 | 2.34 × 10^{−5} | −0.419 707 | 1.87 × 10^{−5} | −0.419 986 | 1.34 × 10^{−4} | −0.436 820 | 1.51 × 10^{−4} | −0.438 127 |

3000 | 1.38 × 10^{−5} | −0.420 246 | 1.10 × 10^{−5} | −0.420 386 | 9.24 × 10^{−5} | −0.439 855 | 9.59 × 10^{−5} | −0.440 324 |

4000 | 1.10 × 10^{−5} | −0.420 480 | 7.78 × 10^{−6} | −0.420 570 | 7.14 × 10^{−5} | −0.441 211 | 6.92 × 10^{−5} | −0.441 458 |

5000 | 9.26 × 10^{−6} | −0.420 609 | 5.65 × 10^{−6} | −0.420 672 | 5.55 × 10^{−5} | −0.441 927 | 5.41 × 10^{−5} | −0.442 146 |

6000 | 6.69 × 10^{−6} | −0.420 686 | 4.44 × 10^{−6} | −0.420 735 | 4.50 × 10^{−5} | −0.442 439 | 4.37 × 10^{−5} | −0.442 607 |

7000 | 5.58 × 10^{−6} | −0.420 745 | 3.44 × 10^{−6} | −0.420 776 | 3.61 × 10^{−5} | −0.442 792 | 3.66 × 10^{−5} | −0.442 933 |

8000 | 4.55 × 10^{−6} | −0.420 774 | 2.69 × 10^{−6} | −0.420 780 | 3.15 × 10^{−5} | −0.443 334 | 3.08 × 10^{−5} | −0.443 173 |

Extrapolated | 0.0 | −0.420 948 | 0.0 | −0.444 784 | ||||

CCSD(T) | −0.398 638 | −0.422 229 | ||||||

CCSDTQ | −0.406 696 | −0.430 244 |

### C. Butadiene

We next consider a larger polyatomic molecule, butadiene. Butadiene has been of continuing interest to quantum chemists, due to the correlation effects in the *π*-conjugated system, and has recently been the target of “exact” methods including *i*-FCIQMC^{88} and composite high-order coupled cluster approaches.^{89}

Unlike in our earlier model arene calculations, here we correlate all electrons except for a frozen 1*s* core. We carry out calculations in the same ANO-L-pVDZ basis^{106} as used in Ref. 88. This generates a space of (22*e*, 82*o*). This is representative of a large “all-electron” exact calculation with the DMRG. We considered both restricted HF canonical orbital as well as MP2 natural orbitals, at the geometry given in Ref. 88. Orbitals were ordered by the GA.

Energies using the different orbital choices up to a modest *M* = 2000 are given in Table V. We observe that the natural orbitals perform better than the Hartree-Fock orbitals, and the split-localized orbitals perform better than the canonical orbitals. Clearly for a molecule of this size locality is already very important. We thus use only the split-localized natural orbitals for the larger *M* calculations on this molecule.

. | Type of orbital . | |||
---|---|---|---|---|

. | . | Canonical . | . | MP2 nature . |

M . | Canonical . | Split-localized . | MP2 natural . | Split-localized . |

500 | −0.500 593 | −0.550 526 | −0.505 492 | −0.552 366 |

1000 | −0.516 424 | −0.554 396 | −0.523 147 | −0.555 348 |

2000 | −0.531 067 | −0.556 172 | −0.537 489 | −0.556 669 |

. | Type of orbital . | |||
---|---|---|---|---|

. | . | Canonical . | . | MP2 nature . |

M . | Canonical . | Split-localized . | MP2 natural . | Split-localized . |

500 | −0.500 593 | −0.550 526 | −0.505 492 | −0.552 366 |

1000 | −0.516 424 | −0.554 396 | −0.523 147 | −0.555 348 |

2000 | −0.531 067 | −0.556 172 | −0.537 489 | −0.556 669 |

The DMRG energies using the split-localized natural orbitals are given in Table VI together with comparison coupled cluster and *i*-FCIQMC energies. We find that we already begin to surpass CCSDT accuracy by about *M* = 2000. By *M* = 6000, we are more than 1 m*E _{h}* below the CCSDT energy, giving an estimate of quadruple and higher order effects. These calculations, however, required considerable resources, especially because there was no point-group symmetry in the split-localized basis: a single

*M*= 3000 sweep in

*C*

_{1}symmetry required 25 h on 42 cores.

M . | Energy . |
---|---|

500 | −0.550 134 |

1000 | −0.554 182 |

2000 | −0.555 899 |

3000 | −0.556 543 |

4000 | −0.556 874 |

5000 | −0.557 050 |

6000 | −0.557 178 |

CCSD(T) | −0.555 002 |

CCSDT | −0.555 959 |

i-FCIQMC | −0.5491(4) |

M . | Energy . |
---|---|

500 | −0.550 134 |

1000 | −0.554 182 |

2000 | −0.555 899 |

3000 | −0.556 543 |

4000 | −0.556 874 |

5000 | −0.557 050 |

6000 | −0.557 178 |

CCSD(T) | −0.555 002 |

CCSDT | −0.555 959 |

i-FCIQMC | −0.5491(4) |

The total energy of butadiene in the ANO-L-pVDZ basis was very recently the subject of a large scale *i*-FCIQMC study. We note that our final DMRG energy is more than 8 m*E _{h}* below the

*i*-FCIQMC energy. Since the DMRG energy is variational, this means that the

*i*-FCIQMC energy is too high, due to initiator bias. The results reported by Daday

*et al.*in Ref. 88 used 1 × 10

^{9}walkers and were limited by the available memory. Clearly, larger numbers of walkers are necessary to remove the initiator bias.

### D. Organometallics

We now turn to some more chemically complex systems, as exemplified by organometallic complexes. As a representative example, we have chosen the oxo-Mn(salen) system (Fig. 10), an analogue of Jacobsen’s catalyst, considered by Ivanic *et al.* using complete active space self-consistent field (CASSCF)^{90,91} and a number of subsequent studies including a recent DMRG-CASSCF study.^{51} We also consider Fe(ii)-porphine, the prototype for biological metalloporphyrins (Fig. 10). Both these molecules are too large to compute all-electron, full basis, DMRG calculations converged to the FCI limit. We thus restrict ourselves here to active-space DMRG calculations. In more realistic studies, such active-space calculations would be augmented by a further more approximate dynamical correlation treatment using perturbation theory,^{8,92,93} configuration interaction,^{55,94} or canonical transformation theory.^{38}

We first consider oxo-Mn(salen). This molecule is used as a simple model for Jacobsen’s epoxidation catalyst.^{95,96} The ordering of the lowest spin states is considered important as different reaction paths have been posited depending on the spin state.^{97} In our DMRG calculations, we used the 6-31G(d) basis set^{81,98,99} and computed the restricted open-shell Hartree-Fock (ROHF) orbitals for the triplet state in the *C*_{1} point group symmetry. Further, the active space was defined by ROHF molecular orbitals with the following dominant atomic orbital characters: (1) five Mn 3*d* orbitals, (2) 2*p _{z}* of the equatorial C, N and O atoms, giving 10

*π*orbitals, (3) 2

*p*and 2

_{x}*p*of the equatorial O and N atoms for the Mn-N and Mn-O

_{y}*σ*bonds, (4) 2

*p*, 2

_{x}*p*and 2

_{y}*p*of the axial Cl atoms, (5) an axial O 2

_{z}*p*, and combination of the axial O 2

_{z}*p*, 2

_{x}*p*orbitals, responsible for the Mn-O

_{y}*σ*and

*π*bonds. The doubly, singly, and unoccupied orbitals were then separately split-localized, and we carried out DMRG(-CI) calculations in this active-space using up to

*M*= 2000 and GA ordered orbitals. (To keep the benchmarks simple and reproducible we do not consider DMRG-CASSCF calculations here.) The results are shown in Table VII. We find that with only 24 orbitals, the active space energy can be converged to tens of microHartrees with a modest

*M*= 2000, a rather modest calculation with this number of orbitals. In this active space, we find that the singlet is lowest in energy, with a very small singlet-triplet energy gap of 0.6 m

*E*. This is reminiscent of the nearly degenerate ground state picture found in previous theoretical works,

_{h}^{51,97}although a quantitative ordering requires augmentation by dynamic correlation.

. | Singlet . | Triplet . | ||
---|---|---|---|---|

M . | DW . | Energy . | DW . | Energy . |

500 | 1.33 × 10^{−5} | −0.304 421 | 1.46 × 10^{−5} | −0.303 777 |

1000 | 4.51 × 10^{−6} | −0.304 648 | 5.20 × 10^{−6} | −0.304 045 |

2000 | 1.19 × 10^{−6} | −0.304 712 | 1.40 × 10^{−6} | −0.304 128 |

. | Singlet . | Triplet . | ||
---|---|---|---|---|

M . | DW . | Energy . | DW . | Energy . |

500 | 1.33 × 10^{−5} | −0.304 421 | 1.46 × 10^{−5} | −0.303 777 |

1000 | 4.51 × 10^{−6} | −0.304 648 | 5.20 × 10^{−6} | −0.304 045 |

2000 | 1.19 × 10^{−6} | −0.304 712 | 1.40 × 10^{−6} | −0.304 128 |

We next discuss the Fe(ii)-porphine system. This is another complex where the correct spin-ordering of the states remains somewhat uncertain.^{92,100,101} Here, we have considered a large active space, consisting of all Fe 3*d* and 4*d* orbitals (the additional 4*d* shell gives the important double-shell correlations) and the full set of ligand *σ* and *π* orbitals. In the *D*_{2h} symmetry, we started from the quintet ROHF orbitals at the triplet geometry described by Groenhof *et al.*^{102} with the cc-pVDZ basis set.^{103–105} We classified orbitals by their dominant atomic or bond character and defined an active space of (44*e*, 44*o*), containing molecular orbitals with the following character: (1) Fe 3*d* and 4*d* orbitals (10 orbitals), (2) 2*p _{z}* orbitals of C and N atoms giving 24

*π*orbitals, (3)

*σ*bonds between Fe 4

*p*, 4

_{x}*p*, and N 2

_{y}*p*and 2

_{x}*p*orbitals (10 orbitals). Finally, all orbitals were split-localized and reordered using the genetic algorithm.

_{y}Table VIII presents the energies obtained using the reverse schedule. The corresponding extrapolations are shown in Fig. 11. We see that although convergence is not rapid in this very large active space, it is still possible to compute all the electronic states to within about 1.5 m*E _{h}* of the estimated exact result at

*M*= 4000. The difficulty of the calculation is similar to that of the all-electron Cr

_{2}calculation, which correlated 48 electrons in 40 orbitals. Within this active space, we find that the energy order is triplet < quintet < singlet. The largest calculations (in

*C*

_{1}symmetry due to split localization) with

*M*= 4000, e.g., for the singlet state, took 15 h per sweep with 40 cores, and are representative of expensive active space DMRG calculations.

. | Singlet . | Triplet . | Quintet . | |||
---|---|---|---|---|---|---|

M . | DW . | E (E_{h})
. | DW . | E (E_{h})
. | DW . | E (E_{h})
. |

1000 | 1.18 × 10^{−4} | −0.833 592 | 1.11 × 10^{−4} | −0.891 395 | 1.05 × 10^{−4} | −0.871 543 |

2000 | 6.42 × 10^{−5} | −0.836 259 | 5.62 × 10^{−5} | −0.893 772 | 4.90 × 10^{−5} | −0.873 551 |

3000 | 4.18 × 10^{−5} | −0.837 218 | 3.61 × 10^{−5} | −0.894 605 | 3.17 × 10^{−5} | −0.874 250 |

4000 | 2.60 × 10^{−5} | −0.837 683 | 2.42 × 10^{−5} | −0.895 003 | 1.95 × 10^{−5} | −0.874 581 |

Extrapolated | −0.839 108 | −0.896 092 | −0.875 330 |

. | Singlet . | Triplet . | Quintet . | |||
---|---|---|---|---|---|---|

M . | DW . | E (E_{h})
. | DW . | E (E_{h})
. | DW . | E (E_{h})
. |

1000 | 1.18 × 10^{−4} | −0.833 592 | 1.11 × 10^{−4} | −0.891 395 | 1.05 × 10^{−4} | −0.871 543 |

2000 | 6.42 × 10^{−5} | −0.836 259 | 5.62 × 10^{−5} | −0.893 772 | 4.90 × 10^{−5} | −0.873 551 |

3000 | 4.18 × 10^{−5} | −0.837 218 | 3.61 × 10^{−5} | −0.894 605 | 3.17 × 10^{−5} | −0.874 250 |

4000 | 2.60 × 10^{−5} | −0.837 683 | 2.42 × 10^{−5} | −0.895 003 | 1.95 × 10^{−5} | −0.874 581 |

Extrapolated | −0.839 108 | −0.896 092 | −0.875 330 |

Many interesting organometallic complexes involve multiple open shell transition metal species. What are the prospects of extending DMRG calculations to such systems? The largest DMRG calculations to date involve 4 open shell transition metal centers with bridging ligands (such as in the Mn_{4}Ca cluster of the oxygen evolving complex in Ref. 52, or the [4Fe-4S] biological redox cofactor in Ref. 53). While such systems are at the frontier of DMRG calculations today, as Ref. 53 shows, calculations with *M* up to 7500, together with extrapolation, yield energies accurate to about 1 m*E _{h}*, within chemical accuracy.

## IV. CONCLUSIONS

In this work, we presented an overview of the theory and practice of the *ab-initio* density matrix renormalization group. In modern implementations, such as in the Block code, the *ab-initio* DMRG is implemented in a black-box manner so that it can be used in a fashion similar to other quantum chemistry methods. The user needs only to specify the desired target number of renormalized states, and the choice of active space and orbitals. With these specifications, all other aspects of the DMRG calculation, such as orbital ordering, perturbative noise, sweep schedule, convergence thresholds, and extrapolations, can be handled automatically by the program.

We have examined the behaviour of the DMRG from the user perspective in a variety of different molecular settings: arenes, diatomics, polyatomics, and organometallics. We summarize our recommendations thus, which are as follows:

For most systems, it is important to balance spatial and energy locality in the orbitals by using split-localized orbitals.

Fiedler vector orbital orderings (default in the Block code) perform well.

For large molecules, we can reason about the number of renormalized states

*M*required to achieve a given accuracy in terms of the effective “width” of the system.For dynamic correlation,

*M*scales roughly linearly with the number of virtuals (basis size) for fixed accuracy.All-electron basis calculations with

*n*of around 2000 (e.g., 20 electrons in 100 orbitals, or 40 electrons in 50 orbitals), converged to chemical accuracy or better and without symmetry, can be considered accessible with cluster computational resources with 50 or so computational cores, with sufficient memory and disk. (Of course, larger calculations are possible with more resources!)_{elec}n_{orb}

We have focused here on relatively simple benchmarks, but of course, many more applications of the DMRG to different systems can be envisaged. We have not explicitly discussed here, for example, excited states, or large polymetallic bioinorganic clusters, although these have been significant targets of study with the DMRG in the literature.^{42,50,52,53} With the simple recommendations above, we believe such systems can now be studied not only by the specialist, but by the general user.

## Acknowledgments

This work was primarily supported by NSF Grant No. OCI-1265278. Additional funding was provided by NSF Grant No. CHE-1265277. The authors acknowledge Helen van Aggelen for useful discussions regarding plot designs, Peter Knowles and Andy May for their help with implementing BLOCK within MOLPRO, Kantharuban Sivalingam for his help with implementing BLOCK within ORCA, and Yihan Shao, for his help with implementing BLOCK within Q-CHEM.

## REFERENCES

*ab initio*programs 2012, see http://www.molpro.net.