This paper describes an interface between the density matrix renormalization group (DMRG) method and the complete active-space self-consistent field (CASSCF) method and its analytical gradient, as well as an extension to the second-order perturbation theory (CASPT2) method. This interfacing allows large active-space multi-reference computations to be easily performed. The interface and its extension are both implemented in terms of reduced density matrices (RDMs) which can be efficiently computed via the DMRG sweep algorithm. We also present benchmark results showing that, in practice, the DMRG-CASSCF calculations scale with active-space size in a polynomial manner in the case of quasi-1*D* systems. Geometry optimization of a binuclear iron-sulfur cluster using the DMRG-CASSCF analytical gradient is demonstrated, indicating that the inclusion of the valence p-orbitals of sulfur and double-shell d-orbitals of iron lead to non-negligible changes in the geometry compared to the results of small active-space calculations. With the exception of the selection of *M* values, many computational settings in these practical DMRG calculations have been tuned and black-boxed in our interface, and so the resulting DMRG-CASSCF and DMRG-CASPT2 calculations are now available to novice users as a common tool to compute strongly correlated electronic wavefunctions.

## I. INTRODUCTION

Strongly correlated many-body problems in quantum chemistry are often a significant challenge to calculations of the electronic structures of complex molecules. This is especially the case when working with transition metal complexes, which are currently of considerable interest and importance in material science and catalysis chemistry. One such many-body problem is the so-called static electron-correlation effect, which results in a multi-reference character in the wavefunction. In traditional quantum chemistry, these many-body problems are solved using multi-reference methods, in which the multi-reference character of a wavefunction is explicitly incorporated. The complete active-space self-consistent field (CASSCF) method^{1,2} is the most common approach amongst the multi-reference methods. As discussed further on, the principal concept of the CASSCF approach is the introduction of an active-space including several orbitals, within which the multi-reference character is fully incorporated. However, the size of the active-space is strictly limited to a maximum of 14–16 orbitals because the complexity of the system increases exponentially with the number of orbitals. For this reason, if the strong electron correlation results from more than 16 orbitals, computations based on the CASSCF method can no longer be performed. Based on these limitations, the density matrix renormalization group (DMRG) method^{3,4} has recently received significant attention in the field of quantum chemistry.^{5–16} However, the DMRG algorithm scales the polynomial complexity with the size of the active-space (AS), even though it works with the full Hilbert space. Therefore, it would be of significant benefit to enlarge the size of the AS that can be applied in the CASSCF method. It would also be useful to extend the DMRG algorithm to work in conjunction with CASSCF-related methods, such as the second-order perturbation theory (CASPT2),^{17,18} *N*-electron valence perturbation theory (NEVPT2),^{19} and multi-reference CI (MRCI) methods.^{20,21} Several illustrative works focusing on this approach have already been reported and several codes have been developed: DMRG-CASSCF,^{22–25} DMRG-CASPT2,^{26} DMRG-NEVPT2,^{27,28} and DMRG-MRCI.^{29} Although these DMRG codes are available as free software, the DMRG method has still not been widely adopted in quantum chemistry. One reason for this is that, in practice, quantum chemistry-DMRG (QC-DMRG) calculations involve numerous computational settings, including the choice of orbitals (such as canonical or localized or other options), orbital ordering, and, of course, the renormalized base size (*M*), in addition to the standard CASSCF calculation parameters. Since a DMRG wavefunction is orbital variant, the appropriate choice of orbitals and the ordering of these orbitals onto the 1D-lattice system are essential to obtaining an accurate result from the DMRG calculations.^{30–32} However, because determining the optimal orbital ordering is an NP-hard problem, this step has historically been performed simply on the basis of the intuition of researchers during pioneering DMRG studies. Recently, we reported a detailed examination of the parameters used in practical DMRG calculations and provided suitable settings for performing a reliable DMRG calculation.^{33} As a consequence, almost all of these settings can potentially be black-boxed so that reliable DMRG calculations can be made without any professional knowledge concerning the DMRG method itself. In addition to this simplification of practical DMRG calculations, predicting stable and/or transition state geometries is a common but crucial component of quantum chemistry research. In most cases, even if we employ an accurate computational method such as the CASPT2 method to evaluate energies and properties, the molecular geometry is usually obtained from either experimental data or from density functional theory (DFT) calculations. This is because the multi-reference methods often require significant computational resources such that it is very difficult to perform a geometry optimization. Even though geometry optimization is usually performed using the DFT method for various types of molecules, it is sometimes problematic in the case of molecules in which a strong static-correlation effect appears. As an example, the potential energy curve of the chromium dimer has long been discussed based on various computational methods, and it has been concluded very recently that only the DMRG-CASPT2 calculation in conjunction with a large active-space can accurately reproduce this curve.^{26} In this regard, it is important to perform geometric optimizations using accurate computational methods to predict the geometry or at least to verify the geometry obtained from DFT calculations. The geometry optimization based on the DMRG-CASSCF method was first reported by Liu and co-workers.^{34} Because the DMRG method is a variational wavefunction ansatz, it fulfills the Hellmann-Feynman theorem^{34,35} and therefore represents a considerable advantage because the analytical energy gradient can be easily derived and efficiently evaluated. Based on this background, we wished to produce a user-friendly DMRG interface to allow the ready performance of DMRG-CASSCF and DMRG-CASPT2 computations, based on an extension of our DMRG code (BLOCK),^{7,15} in which almost all the DMRG settings, except for (*M*), are tuned and black-boxed.

In Sec. II of this report, we provide very brief reviews of the DMRG, CASSCF, and CASPT2 methods and a brief derivation of the analytical energy gradient of a DMRG wavefunction. We also provide a summary explaining that interfacing the DMRG method to the traditional CASSCF/CASPT2 approaches is performed through computations of reduced density matrices (RDMs). Finally in Sec. III, we present some demonstrative applications of the DMRG-CASSCF and DMRG-CASPT2 methods to show their performance with regard to a single-point energy evaluation and also a geometry optimization.

## II. THEORY

In this section, we briefly review the theories and computational methods which are involved in this article. First, we provide a short introduction to matrix product states (MPSs) as a wavefunction representation and show that the DMRG algorithm is derived from the variational principle applied to an MPS wavefunction. Second, we demonstrate that the CASSCF/CASPT2 methods can be linked with the DMRG algorithm through RDMs to perform DMRG-CASSCF/CASPT2 calculations. Third, we note that the energy first derivative of the DMRG is exactly the same as that of the standard CI, which also depends only on RDMs. Finally, we present our strategy for efficiently computing RDMs from the DMRG algorithm.

### A. Matrix product states (MPSs) and the density matrix renormalization group (DMRG) algorithm

Here we present a modern derivation of the DMRG algorithm from the MPS representation of a wavefunction.^{36–38} Although this approach initially seems quite different from the original DMRG formulation using the *block and decimation* scheme,^{5,7} the same formalism eventually results, but from a more natural approach, to quantum mechanics.

Let us start from a full-CI wavefunction representing a linear combination of all possible electron configurations

where $ni={|\u2212\u27e9,|\u2191\u27e9,|\u2193\u27e9,|\u2191\u2193\u27e9}$ and *k* is the number of orbitals being considered. Thus, $|n1\cdots nk\u27e9$ describes a certain electron configuration or, in other words, a Slater determinant. Note that the full-CI wavefunction will consist of an exponentially large number of electron configurations depending on the quantity of orbitals (*k*).

We then decompose the coefficient $Cn1\u2026nk$ into a sequential product of matrices ${\mathbf{A}ini}$. As a result, the wavefunction can be re-written as

where we use the notation $\bm{n}=(n1,\u2026,nk)$ for the sake of simplicity. It should be pointed out that the first and last terms in this expression ($\bm{A}1n1$ and $\bm{A}knk$) are row and column vectors, respectively, and therefore the contraction of all these elements generates a scalar number which reproduces the CI coefficient. To assist the reader in understanding this concept, a graphical representation of an MPS wavefunction is presented in Figure 1. With no approximation, the maximum dimension of these matrices is 4^{k/2}, which is equal to the square root of the Fock-space dimension.

With regard to this MPS representation, it is important that, even if we truncate the size of the matrices, the wavefunction is spanned by the entire collection of electron configurations (that is, the Fock space). In this sense, the MPS mathematically encodes a many-body wavefunction into a very compact formalism. Again for the sake of simplicity, we will assume that each of the matrix objects has the same dimension $M\u07d9\xd7M$ except for the first ($1\xd7M$) and last ($M\xd71$) terms. Consequently, the total degrees of freedom in this wavefunction scale as $O(4M2k)$.

Once we determine a representation for the wavefunction, the ground state can be obtained from a variational principle, applying the parameters of the wavefunction. In this manner, the DMRG algorithm is indeed derived from the variational principle for an MPS wavefunction.

Since it is necessary to compute the overlap between two wavefunctions and the expected value of the Hamiltonian, we first present our approach to efficiently computing these values on the basis of MPS wavefunctions. This occurs in two steps:

- determining the overlap between two MPSs (see also Figure 2)which we can perform with $O(M3k)$ complexity and(3)$\u27e8\Psi A|\Psi B\u27e9=\u2211\bm{n},\bm{n}\u2032(\bm{A}1n1\cdots \bm{A}knk)\u2020\bm{B}1n1\u2032\cdots \bm{B}knk\u2032\u27e8\bm{n}|\bm{n}\u2032\u27e9=\u2211\bm{n}(\bm{A}1n1\cdots \bm{A}knk)\u2020\bm{B}1n1\cdots \bm{B}knk=\u2211\bm{n}\bm{A}knk\u2063\u2020(\cdots (\bm{A}1n1\u2063\u2020\bm{B}1n1)\cdots )\bm{B}knk,$
- determining the expected value of the Hamiltonian (see also Figure 3). As an example, the expected value of an operator acting on the
*i*and*j*orbitals, such that $hija^i\sigma \u2020a^j\sigma $, where*σ*is an index for the electronic spin, is computed aswhere $[a^i\sigma \u2020]nini\u2032=\u27e8ni|a^i\sigma \u2020|ni\u2032\u27e9$ and $[a^j\sigma ]njn\u2032j=\u27e8nj|a^j\sigma |nj\u2032\u27e9$. Since the quantum chemistry Hamiltonian consists of one- and two-body interactions,(4)$\u27e8\Psi 0|hija^i\sigma \u2020a^j\sigma |\Psi 1\u27e9=\u2211\bm{n},\bm{n}\u2032(\bm{A}1n1\cdots \bm{A}knk)\u2020\bm{B}1n1\u2032\cdots \bm{B}knk\u2032hij\u27e8\bm{n}|a^i\sigma \u2020a^j\sigma |\bm{n}\u2032\u27e9=\u2211\bm{n},ni\u2032,nj\u2032hij(\bm{A}1n1\cdots \bm{A}ini\cdots \bm{A}jnj\cdots \bm{A}knk)\u2020\xd7\bm{B}1n1\cdots [a^i\sigma \u2020]nini\u2032\bm{B}ini\u2032\cdots [a^j\sigma ]njnj\u2032\bm{B}jnj\u2032\cdots \bm{B}knk,$we can compute the expected value of the Hamiltonian with naive $O(M3k)\xd7O(k4)$ complexity. It should be noted that optimal computation scaling is achieved by introducing(5)$H^=\u2211ij\sigma hija^i\sigma \u2020a^j\sigma +12\u2211ijkl\sigma \sigma \u2032(ij|kl)a^i\sigma \u2020a^k\sigma \u2032\u2020a^l\sigma \u2032a^j\sigma ,$*complementary operators*to acquire the partial sum of two-body interactions and of $O(M3k3+M2k4)$.^{5,7}

We can then use the Lagrange multiplier technique to find the ground state for which the normalization condition holds

Now, taking into account the variation with respect to the matrix element $\bm{A}ini$ leads to

To assist the reader in understanding this process, $\bm{H}inini\u2032$ and $\bm{S}inini\u2032$ are defined graphically in Figures 4 and 5, respectively. It is of considerable importance that this is an equation for a specific matrix object that is formally independent of the other matrices. Contributions from the other matrices are through the effective Hamiltonian $\bm{H}nini\u2032$ and the overlap matrix $\bm{S}nini\u2032$. Therefore, we can solve a set of equations for each matrix object in an iterative manner until a self-consistent field is obtained. This is done using a *sweep algorithm* in which the effective Hamiltonian (and the overlap matrix) is updated iteratively from left to right and from right to left along the 1D-sequence of the matrix objects.

Since we have only introduced here the basic aspects and most important features of the MPS and DMRG, it is highly recommended that readers refer the original papers^{3,4} and/or several excellent reviews^{14,38–40} for more details regarding MPS and DMRG.

### B. Complete active-space self-consistent field (CASSCF) method

In this subsection, we present a brief overview of the CASSCF method that addresses the manner in which it is connected to the DMRG algorithm, a task that is very easily performed in terms of RDMs.^{41} This is followed by a brief overview of the CASPT2 method, which we will see can also be formulated in terms of RDMs, for which reason we focus on writing formulations of both the CASSCF and the CASPT2 methods in terms of RDMs. It should be noted that this paper uses an orbital notation system which is commonly employed in this field: indices $p,q,r,\u2026$ for any orbital, $i,j,k,\u2026$ for inactive orbitals, $t,u,v,\u2026$ for active orbitals, and $a,b,c,\u2026$ for secondary orbitals.

The CASSCF method was originally proposed by Roos in 1980,^{1} and it has since become established as a standard tool for the investigation of strongly correlated electronic ground and excited states in quantum chemistry. The basic approach associated with the concept of CAS is to define an orbital subspace which can describe all or at least the main part of the strong electron correlation effect, typically referred to as the AS, and to fully incorporate this effect by solving the full-CI problem within this AS. Thus, the *so-called* CAS-CI wavefunction has the form

where $|0\u27e9$ is a reference wavefunction (normally the Hartree-Fock wavefunction) and ${E^I}$ and {*C*_{I}} are all possible excitations and their coefficients within the AS.

Although a full-CI wavefunction is invariant for any orbital rotation, a CAS-CI wavefunction will vary, and so an associated set of orbital rotations must also be optimized. Applying a variational principle to the orbital rotations within the CAS-CI wavefunction leads to the CASSCF method, and the best ground state is obtained by optimizing the CI coefficients and the orbital rotations of the CAS-CI wavefunction. The effects of orbital rotations on the CAS-CI wavefunction are parameterized using a set of single excitations. By introducing a spin-averaged single excitation operator, $E^pq=\u2211\sigma a^p\sigma \u2020a^q\sigma $, the rotated wavefunction can be written as

where

By definition, {*x*_{pq}} is a skew-symmetric matrix, and therefore the term $e\bm{x}$ results in a unitary rotation of the entire orbital space. In practice, such orbital rotations are determined using quadratic convergence algorithms,^{1,42} and a practical CASSCF algorithm therefore alternates CAS-CI calculations with orbital rotations until the energy value converges. Taking the derivative of the expected energy value $\u27e8\Psi 0\u2032|H^|\Psi 0\u2032\u27e9$ with respect to *x*_{pq}, we obtain

where *F*_{pq} is a generalized Fock matrix,

Therefore, the orbital gradient can be evaluated through the 1RDM ($\Gamma qp$) and 2RDM ($\Gamma qspr$), which are defined as

A quadratic convergence algorithm requires at least an approximation of the orbital Hessian matrix and we note here that previous works have concluded that this can also be evaluated using the 1RDM and 2RDM.^{42} At this point, the remaining task is to compute the 1RDM and 2RDM based on a DMRG wavefunction to perform a DMRG-CASSCF calculation.

### C. CASSCF followed by second-order perturbation theory (CASPT2)

Second-order perturbation theory applied to the CASSCF wavefunction (CASPT2)^{17,18} is often employed to incorporate the remaining electron correlation effect outside the AS, which enables one to obtain quantitatively correct results. In this subsection, we will show that the CASPT2 method can also be expressed in terms of RDMs such that the DMRG-CASPT2 method is readily derived. The first-order interacting space of the CASSCF wavefunction is a subspace that consists of single and double excitations from the $|\Psi 0\u27e9$, denoted as *V*_{SD}. According to previously published work, this involves eight orthogonal components,^{17,18}

The first-order wavefunction is expanded by many-body basis functions ${|J\u27e9}\u2208VSD$ as

where ${CJ;J=1\cdots M}$ is a coefficient vector obtained by solving the first-order equation

From the original formulation, the zeroth order Hamiltonian $H^0$ has the form

where $P^0$, $P^K$, $P^SD$, and $P^TQ$ are projectors onto *V*_{0}, *V*_{K}, *V*_{SD}, and *V*_{TQ}, respectively. Note that the one-body operator $F^$ is defined as

where

Consequently, Eq. (17) can be written for each perturbation *V*,

where $SJI=\u27e8J|I\u27e9$ and $VJ=\u27e8J|H^|\Psi 0\u27e9$. Note that the operator $F^$ can be decomposed into diagonal $F^D$ and off-diagonal $F^N$ parts for the sake of convenience

The CASSCF wavefunction satisfies the generalized Brillouin theorem, and therefore $H^0$ must be highly sparse; *fpq* is zero in the case that one index is from the internal subspace and the other is from the secondary subspace. A further simplification can be made by transforming each subspace of *fpq* to a diagonal form, as in

where the ${\u03f5p}$ terms are the eigenvalues for each subspace matrix. This technique represents a *pseudo-canonical* transformation. Because the CASSCF wavefunction is invariant with respect to orbital rotations within each inactive, active, or secondary subspace, it is also invariant with respect to this pseudo-canonical transformation. We can see that $F^N$ has numerous zero blocks and therefore an efficient algorithm to solve the first-order Equation (21) can be derived and implemented on the basis of this partitioning.^{43} The remaining challenge is to efficiently evaluate the matrix elements appearing in Eq. (21), which is formally done by computing RDMs up to fourth order. A bottleneck in this process is, of course, the fourth-order components such that

where $\Gamma uxzwtvyw$ is a component of the 4RDM within the AS. Because a very large quantity of memory is required to explicitly store the 4RDM, the CASSCF wavefunction is first transformed on a pseudo-canonical basis, after which Eq. (24) is evaluated directly, in practical implementations of the CASPT2 method.

### D. DMRG-CASPT2 with the cumulant approximation

In this subsection, we discuss linking the DMRG algorithm with the CASPT2 method. It should be noted that, in the case of a rigorous discussion, the zeroth order Hamiltonian depends on the *M* value of the DMRG wavefunction. However, here we assume that a sufficiently large *M* is employed that the dependency on *M* is converged. A significant issue arises in the practical applications of the DMRG-CASPT2 algorithm, in that the DMRG wavefunction is “variant” with respect to active-active orbital rotation. Therefore, the pseudo-canonical transformation cannot be performed for a DMRG-CASSCF wavefunction without increasing the dimension, *M*, of ${\bm{A}ni}$. There are two possible solutions to this problem: (i) carry out the DMRG calculation once again after the pseudo-canonical transformation or (ii) explicitly compute RDMs and transform these on a pseudo-canonical basis. The first approach leads to the most efficient algorithm for directly computing the components of the first-order Eq. (21). However, the accuracy of the DMRG wavefunction is decreased in this case. That is, a larger *M* would be required to maintain accuracy since, in practice, canonical orbitals are not the best choice for DMRG calculations.^{33} In contrast, the second approach requires using the 4RDM and the transformation of this RDM on a pseudo-canonical orbital basis, which involves the use of a large amount of memory to explicitly store the 4RDM. To avoid this complication, we instead transformed the 1-3 RDMs and employed the cumulant reconstruction method to evaluate the 4-RDM component from the 1-3 RDMs (prefixed by “cu4”).^{26,27,29} This approach should have significant practical benefits, since localized molecular orbitals (LMOs) can be incorporated into the DMRG calculations, and the DMRG algorithm often works very well with such orbitals.^{33} The cumulant reconstruction method for RDMs was originally derived on the basis of spin orbitals, and the previous cu4-DMRG-NEVPT2^{27} and cu4-DMRG-MRCI^{29} methods utilized this spin-orbital approach. In this work, we instead employed the spin-averaged cumulant expansion formalism (see the Appendix),^{44,45} which has advantages for working on the spin-adapted DMRG code. Consequently, the cu4-DMRG-CASPT2 algorithm consists of the following: (i) computing the 1-3 RDMs, (ii) transforming these on a pseudo-canonical basis, and (iii) directly calculating the 4RDM intermediate (Eq. (24) upon reconstruction using Eq. (A1).

### E. Energy first derivative and the Hellmann-Feynman theorem

As noted above, the DMRG algorithm is a variational ansatz for a many-body wavefunction. Therefore, the energy first derivative follows the Hellmann-Feynman theorem for the derivative of each matrix element $\bm{A}ni={A\alpha i\u22121\alpha ini}$. As a result, the energy first derivative of the DMRG wavefunction has exactly the same form as that of the CI wavefunction^{35}

where *a* is an arbitral parameter, such as a Cartesian coordinate, and $hija$ and (*ik*|*jl*)^{a} are the derivatives of the one- and two-particle integrals, respectively,

In this expression, the energy first derivative only depends on the 1RDM and 2RDM, which are obtained from the converged DMRG or the CI wavefunction. One important consequence of the relationship between the DMRG and the CI energy gradients is that the conventional CI gradient code can be utilized to evaluate the DMRG gradient without any modifications.

However, it should be noticed that the coupled-perturbed Hartree-Fock (CPHF) equation must be solved in case that the orbital-SCF is not performed.^{35} In most of the common quantum chemistry packages, the CPHF solver has been implemented for canonical orbitals, while it has not been usually done for localized orbitals which are often employed for the DMRG computations. Thus, the DMRG gradient without the orbital-SCF would be depending on the status of each program package. Furthermore, the DMRG gradient is also influenced by the size of *M* even for the DMRG-CASSCF calculation. Therefore, there is another coupled-perturbed equation for $\bm{A}ini$ which is to be solved by the DMRG response theory.^{46,47} However, this can be ignored if the sufficiently large *M* is used, and in such condition, the DMRG-CASSCF gradient is almost free from the choice of orbitals as we have numerically checked in Sec. III C..

### F. Reduced density matrices from DMRG wavefunctions

To allow for practical implementation, the remaining issue is efficiently computing RDMs from the DMRG wavefunction. Several previous publications have proposed efficient algorithms to compute the 1 and 2RDMs^{41} and the third- or higher-order RDMs,^{28} using the DMRG sweep algorithm.

For the sake of simplicity, we focus here on the 2RDM ($\Gamma jlik$) in the case where $i\u2264j\u2264k\u2264l$, and divide the indices into left (L), dot (D), and right (R) blocks so as to minimize the total computational complexity associated with the DMRG sweep (see Figure 6) as follows:

Here the renormalized operators $[O^]L$ and $[O^]R$ are $M\xd7M$ matrices and $[O^]D$ is a $4\xd74$ matrix. This approach requires $O(M2k2)$ memory to store $[a^k\sigma \u2032\u2020a^l\sigma \u2032]$ and $O(M3k4)$ computational complexity to evaluate all the elements of the 2RDM, just as in the case of the ground state DMRG calculation. Here we do not address the details of the algorithm for the third- and higher-order RDMs, since this would involve a very lengthy discussion. The same concept applied to the 2RDM can also be employed with respect to the higher-order RDMs.^{28}

## III. DEMONSTRATIVE APPLICATIONS

We implemented the interface to our DMRG code in the MOLCAS program package.^{43} This interface connects the DMRG code with the CASSCF, the CASPT2, and the CASSCF analytical gradient codes in the MOLCAS program. Thus, the DMRG-CASSCF calculations either to obtain a single-point energy or for a geometry optimization, as well as the DMRG-CASPT2 calculation of a single-point energy are to be easily performed according to our interface program. In this section, we present some demonstrative applications of the DMRG-CASSCF and cu4-DMRG-CASPT2 methods to molecular systems. First, we consider benchmark calculations for polyene molecules to show the performance and accuracy of the DMRG-CASSCF/cu4-DMRG-CASPT2 methods compared with the standard CASSCF/CASPT2 methods. Second, we further discuss on the accuracy of the cu4-DMRG-CASPT2 method in potential energy curves for the nitrogen dimer as a typical strong electron correlation problem. Finally, we demonstrate the geometry optimization of a [Fe_{2}S_{2}]^{2+} cluster on the basis of the DMRG-CASSCF method.

### A. Polyene molecules

We performed ground state calculations for polyene molecules ranging from butadiene (*n*=4) to eicosadecaene (*n*=20) using the CASSCF/CASPT2 and DMRG-CASSCF/DMRG-CASPT2 methods, with all the *π*-valence orbitals included in the AS. Based on this approach, the number of carbon atoms (*n*) in each molecule was equal to the size of the AS. In addition, for both carbon and hydrogen atoms, Dunning’s correlation-consistent polarized double-zeta basis set (cc-pVDZ) was employed. In the DMRG-CASSCF calculations, the initial orbitals were based on localized molecular orbitals (LMOs) within the AS, where we simply localized all the *π*-valence orbitals together to obtain fully localized orbitals, that is, atomic *p*-like orbitals on each carbon center. Interestingly, the localization of the orbitals was maintained in the converged DMRG-CASSCF wavefunction.

Figure 7(a) shows the computational scaling as a function of the size of the AS for the standard CASSCF and DMRG-CASSCF calculations (with *M* = 200). As we expected, these scaling curves indicate that the DMRG-CASSCF method is advantageous as the size of the AS exceeds 14. The roughly straight line generated in the log-plot of the DMRG-CASSCF data indicates that the DMRG-CASSCF scales polynomially according to the size of the AS in the case of quasi-1*D* systems such as polyene molecules. It should be noted that computational cost of the DMRG method still increases exponentially in the cases of 2*D* and more general systems.^{33} Figure 7(b) and Table I show the energy errors between the DMRG-CASSCF and standard CASSCF results, while Figure 7(c) and Table II show the same data for the cu4-DMRG-CASPT2 and standard CASPT2 calculations. These data are based on using renormalized state sizes (*M*) of 20, 40, 100, 200, and 400. From Figure 7(b) and Table I, it is evident that the DMRG-CASSCF energy converges to the standard CASSCF energy as *M* increases. The energy error achieves chemical accuracy ($<1mEh$) even when *M* is relatively small, since the DMRG works extremely well for polyene molecules. In contrast, it appears that energy errors are more pronounced in the case of CASPT2 calculations. Figure 7(c) and Table II show that the error associated with the calculations using the smallest *M* (*M* = 20) mainly results from the energy error in the DMRG-CASSCF calculations. However, when *M* is larger than 20, the energy error does not depend on either *M* or the size of the AS. There are two possible explanations for this outcome: accurate RDMs can be obtained from the DMRG with relatively small *M* values and/or the loss of accuracy in the DMRG-CASSCF results can be partly recovered by the perturbation. However, it should also be noticed that even with the large *M* value, the cu4-DMRG-CASPT2 results do not match exactly with the CASPT2 results. This is due to inexactness of 4RDM from cumulant approximation, and because the error seems to be worse in the longer polyenes, it indicates that the breakdown of the cumulant approximation will occur for strongly correlated systems. From these results, we conclude that in the case where strong electron correlation is not serious, the DMRG-CASSCF/DMRG-CASPT2 methods are advantageous when computing energy values and also allowing the use of fewer computational resources with minimal loss of accuracy in large AS calculations, while in the case of strongly correlated systems, breakdown of cumulant approximation must be taken care. In Subsection III B, we will further discuss on the performance of the cu4-DMRG-CASPT2 method with one of the strongly correlated systems.

n . | M = 10 . | M = 20 . | M = 40 . | M = 100 . | M = 200 . | M = 400 . |
---|---|---|---|---|---|---|

8 | $3.25\xd710\u22121$ | $2.13\xd710\u22122$ | $4.93\xd710\u22124$ | $2\xd710\u22127$ | $\u22125\xd710\u22127$ | $1\xd710\u22127$ |

10 | $5.54\xd710\u22121$ | $4.94\xd710\u22122$ | $1.29\xd710\u22123$ | $4.7\xd710\u22126$ | $6\xd710\u22127$ | $3\xd710\u22127$ |

12 | $9.67\xd710\u22121$ | $1.06\xd710\u22121$ | $4.87\xd710\u22123$ | $4.89\xd710\u22125$ | $1.1\xd710\u22126$ | $7\xd710\u22127$ |

14 | 1.46 | $1.81\xd710\u22121$ | $1.09\xd710\u22122$ | $1.94\xd710\u22124$ | $5.60\xd710\u22126$ | $1.0\xd710\u22126$ |

n . | M = 10 . | M = 20 . | M = 40 . | M = 100 . | M = 200 . | M = 400 . |
---|---|---|---|---|---|---|

8 | $3.25\xd710\u22121$ | $2.13\xd710\u22122$ | $4.93\xd710\u22124$ | $2\xd710\u22127$ | $\u22125\xd710\u22127$ | $1\xd710\u22127$ |

10 | $5.54\xd710\u22121$ | $4.94\xd710\u22122$ | $1.29\xd710\u22123$ | $4.7\xd710\u22126$ | $6\xd710\u22127$ | $3\xd710\u22127$ |

12 | $9.67\xd710\u22121$ | $1.06\xd710\u22121$ | $4.87\xd710\u22123$ | $4.89\xd710\u22125$ | $1.1\xd710\u22126$ | $7\xd710\u22127$ |

14 | 1.46 | $1.81\xd710\u22121$ | $1.09\xd710\u22122$ | $1.94\xd710\u22124$ | $5.60\xd710\u22126$ | $1.0\xd710\u22126$ |

n . | M = 10 . | M = 20 . | M = 40 . | M = 100 . | M = 200 . | M = 400 . |
---|---|---|---|---|---|---|

8 | $3.58\xd710\u22121$ | $2.96\xd710\u22122$ | $7.29\xd710\u22123$ | $\u22126.07\xd710\u22123$ | $1.23\xd710\u22122$ | $\u22125.04\xd710\u22123$ |

10 | $6.49\xd710\u22121$ | $4.47\xd710\u22122$ | $\u22128.48\xd710\u22124$ | $\u22121.43\xd710\u22122$ | $4.22\xd710\u22123$ | $2.07\xd710\u22123$ |

12 | 1.06 | $1.24\xd710\u22121$ | $2.58\xd710\u22124$ | $\u22122.65\xd710\u22122$ | $\u22122.69\xd710\u22122$ | $\u22122.60\xd710\u22122$ |

14 | 1.59 | $1.96\xd710\u22121$ | $\u22121.50\xd710\u22122$ | $\u22123.75\xd710\u22122$ | $\u22123.66\xd710\u22122$ | $\u22123.62\xd710\u22122$ |

n . | M = 10 . | M = 20 . | M = 40 . | M = 100 . | M = 200 . | M = 400 . |
---|---|---|---|---|---|---|

8 | $3.58\xd710\u22121$ | $2.96\xd710\u22122$ | $7.29\xd710\u22123$ | $\u22126.07\xd710\u22123$ | $1.23\xd710\u22122$ | $\u22125.04\xd710\u22123$ |

10 | $6.49\xd710\u22121$ | $4.47\xd710\u22122$ | $\u22128.48\xd710\u22124$ | $\u22121.43\xd710\u22122$ | $4.22\xd710\u22123$ | $2.07\xd710\u22123$ |

12 | 1.06 | $1.24\xd710\u22121$ | $2.58\xd710\u22124$ | $\u22122.65\xd710\u22122$ | $\u22122.69\xd710\u22122$ | $\u22122.60\xd710\u22122$ |

14 | 1.59 | $1.96\xd710\u22121$ | $\u22121.50\xd710\u22122$ | $\u22123.75\xd710\u22122$ | $\u22123.66\xd710\u22122$ | $\u22123.62\xd710\u22122$ |

### B. Nitrogen dimer

Second, we computed potential energy curves along with triple bond cleavage of the nitrogen dimer on the basis of the cu4-DMRG-CASPT2 method. The dissociation curve of a nitrogen dimer is one of the typical strong electron correlation problems in quantum chemistry, and therefore, it is of great use to examine methods for quantum many-body problems. Here, we would like to discuss whether or not cumulant approximation works and mention a guideline for practical DMRG-CASPT2 applications. Figure 8 shows energy errors in the DMRG-CASSCF and cu4-DMRG-CASPT2 calculations from the standard CASSCF and CASPT2 results for three different bond lengths of the nitrogen dimer (*R*_{N–N}).

In the DMRG-CASSCF calculations, the energy well converges to the standard CASSCF energy as *M* increases. While in the cu4-DMRG-CASPT2 calculations, there are still non-negligible energy differences from the standard CASPT2 results even though the energy converges with respect to *M*. As a consequence, the remaining error can be understood as an error from the cumulant approximation. Although the error is up to 1 m*E*_{h} in this case, it seems to be considerable in the intermediate bond length, indicating that the error from the cumulant approximation is getting worse for a strongly correlated system.

It should be noted that cumulant-approximation-free dynamical correlation methods based on the DMRG-CASSCF wavefunction were also proposed previously.^{48–50} In all these methods, a key to perform practical multi-reference calculations coupled with the DMRG-CASSCF wavefunction is to avoid the explicit computation of the 4RDM. However, as we have also mentioned in Sec. II D, this is not trivial with the LMO basis that is often employed in the DMRG part of computation. Therefore, the cumulant-based approximation is still of great use in practical DMRG-CASPT2 calculations although we must take care of the expected breakdown in cases of extremely correlated systems.

### C. Iron-sulfur cluster

Finally, we present an illustrative application to the determination of the energy first derivative of a DMRG-CASSCF wavefunction. For this example, we employ a simple iron-sulfur cluster [Fe_{2}S_{2}]^{2+} and perform a geometry optimization. In the case of *D*_{2h} symmetry, there are only two parameters (that is, the Fe–Fe and Fe–S distances) and these are independent, although it is still helpful to examine whether the strong electron correlation affects the geometry of this cluster. Table III summarizes the optimized geometries based on DFT and DMRG-CASSCF calculations. In the DFT geometry optimizations, we examined two common functionals, B3LYP and BP86, which have often been employed in computational studies of iron complexes and clusters. In the DMRG-CASSCF geometry optimizations, we also assessed several AS choices: (1) an AS including 3d orbitals from the two iron centers, resulting in 10 electrons and 10 orbitals, denoted as (10e, 10o), (2) including 4d orbitals in addition to (10e, 10o), resulting in a (10e, 20o) AS, (3) including the 3p orbitals of the two sulfur atoms in addition to the (10e, 10o), resulting in a (22e, 16o) AS, and (4) including 4d and 4s orbitals in addition to the (22e, 16o), resulting in a (22e, 28o) AS. The optimized geometries obtained from the two functionals were computed and were found to be very similar to one another in the high-spin state, although the two approaches gave somewhat different Fe–Fe distances in the low-spin state. The energy gaps between the high-spin and low-spin states (the H–L gaps) were also considerably different between these two functionals. Thus, it is concluded that the effect of the functional is quite large, at least in terms of energy, for a metal cluster such as this. The optimized geometry from the standard CASSCF (10e, 10o) calculation gave a result similar to that obtained from the B3LYP functional, except for the Fe–Fe distance in the high-spin state. The DMRG-CASSCF (10e, 20o) result was also very close to that of the CASSCF (10e, 10o), even though the low-spin state was not converged to the correct ground state. Interestingly, the DMRG-CASSCF (22e, 16o) calculation gave completely different geometries for the high-spin and low-spin states. In addition, both the Fe–Fe and Fe–S distances were found to be longer compared to the values obtained from the DFT and standard CASSCF (10e, 10o) calculations. These data suggest that the effect of the size of the AS is not converged at the standard CASSCF level of theory. Finally, the DMRG-CASSCF (22e, 28o) calculations, representing the largest AS calculation in this work, gave a geometry in which the Fe–Fe and Fe–S distances were surprisingly similar to those obtained from the DFT calculations. Furthermore, the resulting H–L gap was also different from those generated by the CASSCF and DMRG-CASSCF calculations working with a smaller AS. As a result, the geometry and the H–L gap obtained from the DMRG-CASSCF (22e, 28o) calculations were extremely close to those resulting from the B3LYP calculations. These results demonstrate that geometry optimization via the CASSCF method is also significantly affected by the lack of a large static electron correlation effect in such molecules. Though it is probable that the B3LYP calculations contained a fortunate error cancellation, it is of considerable importance that we were able to guarantee validity of the B3LYP geometry optimization and energy value for this molecule.

. | High-spin . | Low-spin . | . | ||
---|---|---|---|---|---|

. | Fe–Fe (Å) . | Fe–S (Å) . | Fe–Fe (Å) . | Fe–S (Å) . | H–L gap (eV) . |

BP86 | 2.939 | 2.273 | 2.635 | 2.178 | 1.83 |

B3LYP | 2.953 | 2.264 | 2.794 | 2.203 | 0.66 |

CASSCF(10e, 10o) | 2.835 | 2.245 | 2.786 | 2.228 | 0.21 |

DMRG-CASSCF(10e, 20o) | 2.825 | 2.244 | n.c. | n.c. | N/A |

DMRG-CASSCF(22e, 16o) | 3.204 | 2.453 | 3.145 | 2.449 | 0.23 |

DMRG-CASSCF(22e, 28o) | 2.982 | 2.356 | 2.681 | 2.219 | 0.83 |

. | High-spin . | Low-spin . | . | ||
---|---|---|---|---|---|

. | Fe–Fe (Å) . | Fe–S (Å) . | Fe–Fe (Å) . | Fe–S (Å) . | H–L gap (eV) . |

BP86 | 2.939 | 2.273 | 2.635 | 2.178 | 1.83 |

B3LYP | 2.953 | 2.264 | 2.794 | 2.203 | 0.66 |

CASSCF(10e, 10o) | 2.835 | 2.245 | 2.786 | 2.228 | 0.21 |

DMRG-CASSCF(10e, 20o) | 2.825 | 2.244 | n.c. | n.c. | N/A |

DMRG-CASSCF(22e, 16o) | 3.204 | 2.453 | 3.145 | 2.449 | 0.23 |

DMRG-CASSCF(22e, 28o) | 2.982 | 2.356 | 2.681 | 2.219 | 0.83 |

We further discuss dependencies to a choice of orbitals and to *M* value in the geometry optimization. Since orbital gradients inside the AS must be approximately zero in the DMRG-CASSCF wavefunction, gradients should be less influenced by a choice of orbitals. Table IV summarizes the energy and absolute forces acting on Fe and S atoms of [Fe_{2}S_{2}]^{2+} at the BP86-optimized geometry. As *M* increases, either using localized natural orbitals (LNOs) or using (decalized) natural orbitals (NOs), it converges to the same number. Thus, it is numerically elucidated that the gradient from DMRG-CASSCF has almost no dependency to a choice of orbitals in geometry optimization. It should be also noted that gradients also depend on the *M* value, which seems to be related to the convergence in energy. As a result, using LNOs gives the faster convergence both in energy and gradient value. Table V shows dependency to *M* value in two different AS: (22e in 16o) and (22e in 28o). Averaged number of DMRG-CASSCF iterations and approximate computational cost per optimization cycle are also displayed. In principle, when the size of AS is enlarged, it is necessary to increase *M* value. Of importance here is that this is true for the convergence of gradient value as well as the convergence of energy. Approximate computational cost increases moderately by increasing *M* and the size of AS. Therefore, CASSCF geometry optimization with a large AS is of great use in practical calculations for strongly correlated systems.

. | LNOs (C_{s} symmetry)
. | NOs (D_{2h} symmetry)
. | ||||
---|---|---|---|---|---|---|

M
. | E/E_{h}
. | $|fFe|$ . | $|fS|$ . | E/E_{h}
. | $|fFe|$ . | $|fS|$ . |

500 | −3339.366 075 | 0.040 107 | 0.028 364 | −3339.358 088 | 0.041 063 | 0.031 440 |

1000 | −3339.366 151 | 0.039 948 | 0.028 252 | −3339.364 912 | 0.040 125 | 0.028 703 |

2000 | −3339.366 156 | 0.039 937 | 0.028 245 | −3339.366 114 | 0.039 954 | 0.028 265 |

. | LNOs (C_{s} symmetry)
. | NOs (D_{2h} symmetry)
. | ||||
---|---|---|---|---|---|---|

M
. | E/E_{h}
. | $|fFe|$ . | $|fS|$ . | E/E_{h}
. | $|fFe|$ . | $|fS|$ . |

500 | −3339.366 075 | 0.040 107 | 0.028 364 | −3339.358 088 | 0.041 063 | 0.031 440 |

1000 | −3339.366 151 | 0.039 948 | 0.028 252 | −3339.364 912 | 0.040 125 | 0.028 703 |

2000 | −3339.366 156 | 0.039 937 | 0.028 245 | −3339.366 114 | 0.039 954 | 0.028 265 |

22e in 16o . | |||||
---|---|---|---|---|---|

M
. | E/E_{h}
. | Fe–Fe (Å) . | Fe–S (Å) . | CASSCF iters. . | CPU/cycle . |

1000 | −3339.390 763 | 3.145 | 2.449 | 12 | 2 h |

2000 | −3339.390 782 | 3.145 | 2.449 | 5 | 3.5 h |

22e in 28o | |||||

M | E/E_{h} | Fe–Fe (Å) | Fe–S (Å) | CASSCF iters. | CPU/cycle |

1000 | −3339.568 975 | 2.691 | 2.223 | 20 | 5 h |

2000 | −3339.574 003 | 2.681 | 2.219 | 15 | 13 h |

22e in 16o . | |||||
---|---|---|---|---|---|

M
. | E/E_{h}
. | Fe–Fe (Å) . | Fe–S (Å) . | CASSCF iters. . | CPU/cycle . |

1000 | −3339.390 763 | 3.145 | 2.449 | 12 | 2 h |

2000 | −3339.390 782 | 3.145 | 2.449 | 5 | 3.5 h |

22e in 28o | |||||

M | E/E_{h} | Fe–Fe (Å) | Fe–S (Å) | CASSCF iters. | CPU/cycle |

1000 | −3339.568 975 | 2.691 | 2.223 | 20 | 5 h |

2000 | −3339.574 003 | 2.681 | 2.219 | 15 | 13 h |

## IV. CONCLUSIONS

In this article, we briefly reviewed the DMRG algorithm and its connection to the CASSCF and CASPT2 methods and to the energy first derivative of a DMRG wavefunction. It has been shown that the DMRG algorithm is linked to the CASSCF and CASPT2 methods through reduced density matrices. Providing an interface between these methods is highly beneficial, as documented herein using the DMRG interface implemented in the MOLCAS program package. This DMRG interface enables one to easily perform single-point energy calculations by employing the DMRG-CASSCF and DMRG-CASPT2 methods. DMRG-CASSCF can also be employed to perform geometry optimizations.

The results presented herein demonstrate that the computational costs of DMRG-CASSCF calculations scale according to the polynomial complexity, and that improved scaling occurs when the size of the active-space is greater than 14. We have also confirmed that the DMRG-CASSCF energy results converge to those generated using the CASSCF method when *M* is sufficiently large. Of more interest is that the energy errors resulting from DMRG-CASPT2 calculations exhibit a somewhat complicated dependence on *M*. It is also particularly noteworthy that the DMRG-CASPT2 results agree very well with the standard CASPT2 results except at the smallest *M* value, implying that perturbation treatment can recover a portion of the discarded interactions in the DMRG-CASSCF wavefunction. Geometry optimization of the [Fe_{2}S_{2}]^{2+} cluster was performed using the DMRG-CASSCF method to investigate the extent to which the static electron correlation affects the geometry of transition metal complexes. In this process, we examined several choices of the active-space, including up to 28 orbitals, a value which is not possible when using the standard CASSCF approach. It was determined that the geometry is significantly affected by the active-space selection, and that calculations using large active-spaces are required to obtain results that converge. In this demonstration, the B3LYP result agreed very well with the large active-space DMRG-CASSCF result, although there are many challenges associated with using the B3LYP functional for such a complicated molecule. These results support one empirical conclusion, that DFT-optimized geometries are useful even in the case of transition metal complexes if the appropriate functional is selected, although it should not be assumed that the DFT result is always correct. DMRG-based computations can be readily performed using our DMRG interface, and therefore this interface should highly be useful in the computations of accurate geometries and energies for transition metal complexes and clusters which include large active-space static electron correlation problems.

## ACKNOWLEDGMENTS

N.N. would like to thank Professor G. Chan and Professor J. Hasegawa for fruitful discussions and comments. This work was financially supported by JSPS Kakenhi (Grant Nos. JP15K20832 and JP15H03770), by JST-CREST, by Building of Consortia for the Development of Human Resources in Science and Technology, MEXT, Japan, and by the NSF (Grant No. NSF-CHE 1657286).

### APPENDIX: SPIN-AVERAGED CUMULANT EXPANSION OF 4RDM

As an example, the fourth-order cumulant $\Lambda ijklabcd$ (when incorporating a small correction from Ref. 45) has the form

where *i*, *j*, *k*, *l*, *a*, *b*, *c*, *d* are the indices of spatial orbitals, and the numbers in parentheses denote the number of permutations for each term. In this work, the 4RDM was reconstructed from 1 to 3 RDMs as an approximation to $\Lambda ijklabcd\u22480$.