Projection-Based Density Matrix Renormalization Group in Density Functional Theory Embedding

The density matrix renormalization group (DMRG) method has already proved itself as a very efficient and accurate computational method, which can treat large active spaces and capture the major part of strong correlation. Its application on larger molecules is, however, limited by its own computational scaling as well as demands of methods for treatment of the missing dynamical electron correlation. In this work, we present the first step in the direction of combining DMRG with density functional theory (DFT), one of the most employed quantum chemical methods with favorable scaling, by means of the projection-based wave function (WF)-in-DFT embedding. On two proof-of-concept but important molecular examples, we demonstrate that the developed DMRG-in-DFT approach provides a very accurate description of molecules with a strongly correlated fragment.

embedding.On the two proof-of-concept but important molecular examples, we demonstrate that the developed DMRG-in-DFT approach provides a very accurate description of molecules with a strongly correlated fragment.
Strong correlation plays a crucial role in many aspects of chemistry, such as bond breaking processes, open-shell systems, excited electronic states, as well as in catalysis. 1,2Accurate and efficient description of strongly correlated molecules, however, belongs to long-standing challenges of quantum chemistry.In principle, it can be accounted for by the exact full configuration interaction (FCI) method, but it is prohibitively expensive due to its exponential scaling.In order to bypass the limitations of FCI, several approximate polynomially scaling wave function (WF) methods were developed over the years, which can be systematically improved towards FCI.In case of molecules with weakly correlated electrons, such as organic molecules composed from the main elements and at equilibrium geometries, the most prominent example is undoubtedly the coupled cluster method, 3 whereas the concept of the complete active space (CAS) 4 can be considered as a standard tool for strongly correlated molecules, such as transition metal complexes and bond breaking processes.The last two cases are also the focus of this work.
The complete active space self-consistent field (CASSCF) method, 5 which couples FCI in a small active space with orbital optimization, is usually the starting point of multireference (MR) calculations.The missing dynamical electron correlation is then taken into account by post-SCF methods, such as the complete active space second-order perturbation theory (CASPT2), 6 the second-order n-electron valence state perturbation theory (NEVPT2), 7 or the multireference configuration interaction (MRCI). 2 The common hurdle of all these methods is the limited CAS size to less than 20 orbitals, due to the FCI exponential scaling.
Since many molecules, such as transition metal complexes, require larger CAS than FCI can handle, several approximate FCI solvers have been developed, one of them being the density matrix renormalization group (DMRG) method. 81][12] This sparked interest in development of many post-DMRG methods for treatment of the missing (out-of-CAS) dynamical correlation are available. 13However, these WF-based methods are still too costly for large systems of particular interest.Their alternative, the density functional theory (DFT) represent a cost-effective approach applicable to very large molecules, which however, has its own limitations.The major shortcomings of DFT are undoubtedly the approximate form of the exchange-correlation functional as well as the single reference character, which makes it unsuitable for strongly correlated problems. 14e way of extending the range of applicability of accurate (single or multireference) WF-based methods can be achieved by means of the quantum embedding. 15This approach relies on locality of chemical interactions and splits the whole system into the active subsystem that is treated at a high level, and the environment subsystem that is treated at a lower level of theory. 15,16Previously, Neugebauer, Reiher, and co-workers presented the first and to the best of our knowledge the only attempt to embed DMRG calculations in DFT environment by means of the frozen density embedding approach 17 for treatment of strongly correlated systems.However, due to the approximate form of the non-additive kinetic potential (NAKP), their proof-of-principle applications were restricted to systems in which the active subsystem is not covalently bonded to the environment.
The projection-based DFT (PB-DFT) embedding 18 method is free of the NAKP problem, due to the orthogonality of occupied orbitals of both subsystems, which is achieved by the level shift projection operator. 18This additionally ensures that the sum of energies of the active system and the environment effects is equal to the energy of the full system if both fragments are treated at the same level of theory.Encouraged by an impressive performance of the projection-based embedding for various chemical systems such as, transition metal catalysis, enzyme reactivity, or battery electrolyte decomposition, 19,20 as well as by robustness of the DMRG method, herein we develop and implement the DMRG-in-DFT projection-based embedding method.As demonstrated in the remainder of this letter, this approach has a tremendous potential for applications to large strongly correlated systems.
The DMRG method is a variational procedure for approximating the exact FCI wave function with the so called matrix product state (MPS). 21The FCI wave function in the occupation basis representation reads as where occupation of each orbital corresponds to α i ∈ {|0 , | ↓ , | ↑ , | ↓↑ } and the expansion coefficients c α 1 ...αn form the FCI tensor.By successive applications of the singular value decomposition (SVD), the FCI tensor can be factorized to the MPS form 21 where A[j] α j are the MPS matrices specific to each orbital and the newly introduced auxiliary indices i j are contracted over.If the MPS factorization is exact, the dimensions of the MPS matrices grow in a similar fashion as the size of the original FCI tensor, i.e. exponentially (with an increasing system size).In DMRG, the dimensions of auxiliary indices are bounded.
These dimensions are called bond dimensions and are usually denoted with M .
A practical version of DMRG is the two-site algorithm, which provides the wave function in the two-site MPS form For a given pair of adjacent indices [i, (i + 1)], W is a four-index tensor, which corresponds to the eigenfunction of the second-quantized electronic Hamiltonian expanded in the tensor product space of four tensor spaces.The tensor spaces are defined on an ordered orbital chain, so called left block (M l dimensional tensor space), left site (four dimensional tensor space of i th orbital), right site (four dimensional tensor space of (i + 1) th orbital), and right block (M r dimensional tensor space).In Eq. 4, h pq and pq|rs denote standard one and two-electron integrals in the molecular orbital basis, and σ and σ denote spin.The MPS matrices A are obtained by successive application of SVD with truncation on W's and iterative optimization by going through the ordered orbital chain from left to right and then sweeping back and forth. 11The maximum bond dimension (M max ) which is required for a given accuracy, can be regarded as a function of the level of entanglement in the studied system. 22 the following, we will briefly describe the projection-based embedding WF-in-DFT technique.The WF-in-DFT embedding procedure starts with an initial DFT calculation of the whole system.Based on some criteria for associating the molecular orbitals to the active and environment subsystems, the corresponding density matrix γ is partitioned into the active subsystem A and the environment subsystem B, γ A and γ B , respectively.Originally, this was achieved by means of the occupied orbitals localization and Mulliken population analysis, 18 though alternative more robust approaches have also been developed. 23,24In case of the DFT-in-DFT embedded calculation, the total energy can be expressed as 19 where E DFT denotes the DFT energy evaluated using the bracketed density matrix, γ A emb is the embedded subsystem A density matrix, and P B is a projection operator enforcing mutual orthogonalization, P B = Sγ B S. S denotes the atomic orbital overlap matrix.In the limit where the level shift parameter µ → ∞, the A and B orbitals are exactly orthogonal, but µ is for practical purposes taken to be 10 6 , causing negligible error. 18The embedding potential v emb contains all interactions between subsystems A and B The matrix g groups all the two-electron contributions (Coulomb, exchange, and exchangecorrelation). Because, the projection-based embedding approach is free from non-additive kinetic energy problem 18 it is formally exact, i.e. when the active part was treated with the same exchange-correlation functional as the environment, it would be equivalent to the Kohn-Sham solution of the entire system.
The Fock matrix of subsystem A embedded in B has the following form 19 where h is the core Hamiltonian matrix and it is self-consistently optimized with respect to γ A emb .In case of single reference WF-in-DFT calculations, HF-in-DFT with the following effective core Hamiltonian precedes the WF calculation.For MR problems, CASSCF-in-DFT can be performed. 25wever, since we employ the accurate DMRG which approaches the FCI solution of the active subsystem, we are free to use HF-in-DFT for the MR problems.
Most importantly, the DFT-in-DFT method can be straightforwardly employed for a WF-in-DFT embedding where the active subsystem is treated with the DMRG method and the environment subsystem is described with the DFT method.Then the DMRG-in-DFT energy is simply obtained by substituting the DFT energy of the active subsystem A with the DMRG energy as In this equation, E DMRG [Ψ A MPS ] is the DMRG energy of the active subsystem corresponding to the MPS wave function |Ψ A MPS , which minimizes the active subsystem Hamiltonian (4) with the one-electron part replaced by the effective core Hamiltonian from Eq. 8.
The WF-in-DFT embedding method has been implemented in Psi4NumPy quantum chemistry software 26 which was interfaced with the MOLMPS 27 DMRG code.The developed method was then used to study two benchmark problems (see Figure 1) which have a strongly correlated active part coupled to the environment, namely the triple bond stretching in propionitrile (CH 3 CH 2 CN) and the conformational isomerization of the model iron- 28 which is a prototype of a transition metal complex with the non-innocent nitrosyl ligand relevant to medicinal applications. 29Regarding the low-level method, all the DFT calculations employed the B3LYP 30,31 density functional.On the other hand, all the high-level DMRG calculations were warmed-up with the CI-DEAS procedure 11,22 and took advantage of the dynamical block state selection (DBSS), 32 which adjusts the actual bond dimensions to fit the desired (pre-set) truncation error (TRE).The initial DMRG orbital orderings were optimized with the Fiedler method. 33The complementary calculations listed below were carried out in the following programs: CCSD in Psi4, 26 CASSF/DMRG-SCF in Orca, 34 adiabatic connection (AC) in GammCor, 35 and internally contracted MRCI in MOLPRO. 36 our first example, we study the triple bond stretching in propionitrile (CH 3 CH 2 CN) molecule.The equilibrium geometry of propionitrile employed in this work is given in the Supporting Information (SI, Table S1).For the WF-in-DFT calculations, we have employed the cc-pVDZ 37 basis set.The active subsystem comprised the -CN group and the orbitals  were partitioned into both subsystems by means of the SPADE procedure. 23The stretching of the CN bond was probed by the accurate DMRG-in-B3LYP calculations with TRE = 10 −6 .
For comparison, we also carried out the CCSD-in-B3LYP, as well as the CCSD and DMRG calculations for the entire molecule.The frozen-core approximation was employed for the aforementioned DMRG calculations leading to the FCI space of 22 electrons in 77 orbitals and TRE was pre-set to 10 −5 .S2.
As it is known, the CCSD method notoriously fails in describing correctly the triple bond breaking due to its single-determinant nature.It e.g.predicts a nonphysical bump on PES of N 2 molecule in the intermediate stretching region (around 2.2 Å). 38   see in Figure 2, that the situation is unsurprisingly very similar for the triple C-N bond stretching in CH 3 CH 2 CN.The CCSD method provides much higher dissociation energies for the intermediate stretching region, than the frozen-core DMRG (at 2.5 Å, the error is ∼ 2.4 eV).CCSD-in-B3LYP behaves even slightly worse than CCSD itself.On the other hand, there is a huge improvement between CCSD-in-B3LYP and DMRG-in-B3LYP in description of the triple C-N bond stretching process.At 2.5 Å, the error of DMRG-in-B3LYP with respect to DMRG is 0.9 eV, whereas for the CCSD-in-B3LYP method this error is 3.3 eV.The DMRG method as a genuine MR method is able to properly describe this process.The difference between DMRG-in-B3LYP and DMRG, which is essentially very similar to the difference between CCSD-in-B3LYP and CCSD, thus can be attributed to the lower-level (B3LYP) description of the remaining electrons plus errors of the PB-DFT embedding (density-driven errors or errors coming from the non-additivity of the exchangecorrelation energy 39 ).
As our second example, we have studied the conformational isomerization of the model iron-nitrosyl complex [Fe(CN) 5 (NO)] 2− .The B3LYP optimized geometries of the standard, flat, and reversed isomers of [Fe(CN) 5 (NO)] 2− (see Figure 1b) were taken from Ref. 28 (also given in Table S3-S5).For computational reasons, we used the smaller 6-31G 40,41 basis.The active subsystem was formed by [Fe-NO] 3+ and partitioning of the orbitals into subsystems was carried out by means of the SPADE procedure. 23In order to decrease the size of the virtual space, we employed the two-shell concentric localization 42   (16,16).All CASSCF natural orbitals are shown in Figures S3-S9).In the smaller CAS (14,15), we performed CASSCF computations, which were then corrected for the dynamical electron correlation by means of strongly contracted NEVPT2, the adiabatic connection (AC), 43,44 and the linearized-ACintegrand approximation AC0. 43,44The later two have the advantage of favourable scaling with respect to the CAS size and thus represent an ideal choice for approximate FCI solvers such as DMRG. 45In CAS (16,16), we performed the DMRG-SCF calculations with fixed bond dimensions equal to 2000 and subsequent AC/AC0 in order to probe the effect of the missing dynamical electron correlation.
Table 1 shows the natural orbital occupation numbers (NOONs) of the four orbitals around the Fermi level for the largest active space employed, i.e.CAS (16,16) (all occupation numbers can be found in SI).
The occupation numbers largely deviate from 2 (and 0) and confirm the non-innocent nature of the nitrosyl ligand, indicating the significant multireference character of the investigated systems.Moreover, looking at the four aforementioned orbitals (Figures S7-S9), a ∆E S→F denotes the energy difference between flat (F) and standard (S) isomers.b ∆E S→R denotes the energy difference between reverse (R) and standard (S) isomers.
Table 2 shows the reaction energies of three stable isomers involved in the [Fe(CN) 5 (NO)] 2− complex conformational isomerization computed by various single and multi-reference methods as well as with the CCSD and DMRG methods embedded in the HF or DFT environment.The graphical summary is depicted in Figure 3.Because of the significant multireference character in all three isomers, Figure 3 and Table 2 indicate that the single reference methods (B3LYP and CCSD), in contrast to all state-of-the-art multireference approaches, incorrectly predict the reverse isomer to have the highest energy.At the CAS(14,15) level, we can observe that adding the dynamical electron correlation on top of CASSCF by means of NEVPT2 and AC0/AC results in a larger ∆E S→F by 0.6 -0.7 eV, whereas ∆E S→R is affected only slightly.
More importantly, AC0 provides very similar energy gaps as NEVPT2 (within 0.16 eV in case of ∆E S→R ), as was already pointed out previously. 44The canonical AC method captures even more correlation energy than its linearized AC0 approximation and the AC (16,16)   results together with the icMRCISD(4,4) results represent our best estimates of the energy gaps, in particular 1.9-2.14eV for ∆E S→F and ∼1.40 eV for ∆E S→R .
Looking at the results of the embedded calculations in Table 2, one can see that CCSDin-HF as well as CCSD-in-B3LYP underestimate the ∆E S→F gap even more than CCSD and predict incorrectly that the flat isomer is lower in energy than the reverse one (by 0.8 eV and 0.6 eV, respectively).On contrary, the results of the DMRG embedded calculations are in a very good agreement with our best estimates of the energy gaps.Both DMRG-in-HF as well as DMRG-in-B3LYP provide ∆E S→F gaps within the margins of the MR methods, the DMRG-in-B3LYP ∆E S→R gap is slightly lower (by ∼0.2 eV).The DMRG-in-HF method achieves a perfect agreement of both energy gaps with our best estimates obtained by the state-of-the-art MR methods, which confirms that the Fe-NO moiety is mainly responsible for the electronic structure properties of the [Fe(CN) 5 (NO)] 2− complex.
In this letter, we present the projection-based DMRG-in-DFT embedding method and we test its performance on two benchmark problems, namely the triple bond stretching in proach is the size of the virtual space which, even when it is truncated, 42 might be too large for DMRG.It is also the reason, why we were limited to smaller basis sets.However, in case of larger basis sets, the concept of CAS can be used in which DMRG is combined with some post-DMRG method 13 which will the subject of our following works.

Supporting Information
Equilibrium geometry of CH

Figure 2 ,
Figure 2, shows the potential energy surfaces (PES) [differences with respect to minima: E(r CN ) − E min ] corresponding to the triple C-N bond stretching in propionitrile.The results obtained by B3LYP, CCSD, CCSD-in-B3LYP, and DMRG-in-B3LYP are compared against the exact curve obtained by the frozen-core DMRG method.The individual absolute energies are provided in TableS2.

Figure 2 :
Figure 2: Comparison of the individual dissociation energy curves corresponding to the triple C-N bond stretching in CH 3 CH 2 CN.All calculations employ the cc-pVDZ basis set.

CH 3 CH 2
CN and conformational isomerization of [Fe(CN) 5 (NO)] 2− , a prototype of the transition metal complex containing a non-innocent ligand.Both of these systems exhibit a significant multireference character.Our numerical results indicate that the DMRG-in-DFT provides a viable way toward accurate description of molecules containing strongly correlated fragment.In case of the triple bond stretching in CH 3 CH 2 CN, the DMRG-in-B3LYP method substantially outperformed the single-reference CCSD and CCSD-in-B3LYP methods, whereas in case of the [Fe(CN) 5 (NO)] 2− complex, the DMRG-in-B3LYP and DMRGin-HF methods provided the energy gaps between individual isomers that are in very good agreement with the state-of-the-art multireference approaches.This work represents the first step toward combining DMRG with PB-DFT embedding.The biggest bottleneck of this ap- leading to the active subsystem FCI space comprising 38 electrons in 102 orbitals.For comparison, we also carried out the B3LYP and CCSD calculations as well as calculations with different CAS-based MR methods.The smallest CAS(4,4)comprising the two NO π * orbitals together with the Fe 3d xz and 3d yz was employed for internally contracted MRCI with singles and doubles (icMRCISD) calculations.The larger CAS(14,15)contained the NO π (two), π * (two), σ, σ * , and Fe 3d (five), 4d (3 counterparts to the occupied 3d orbitals: 4d xy , 4d xz , and 4d yz ), plus one equatorial σ orbital with the Fe 3d x 2 −y 2 and C 2p x/y contributions.This CAS(14,15)was augmented with one occupied axial orbital of σ character to form CAS

Table 1 :
DMRG-SCF(16,16)Natural Orbital Occupation Numbers for the Individual [Fe(CN) 5 (NO)] 2− Standard (S), Flat (F), and Reverse (R) Isomers.thattheirelectron density is mainly localized to the Fe-NO region, which corroborates the use of the WF-in-DFT embedding, in which the WF method, however, should be able to correctly describe the MR character of the Fe-NO moiety.The strongest MR character is observed for the reverse isomer.In this case, the weight of the HF reference in the DMRG-SCF(16,16)wave function is only 64% and one can expect that the conventional single reference approaches might be inappropriate.

Table 2 :
Reaction Energies in eV Corresponding to the Conformational Isomerization of [Fe(CN) 5 (NO)] 2− Complex Calculated with Different Methods and 6-31G Basis Set.

Table S2 :
Absolute energies of CH 3 CH 2 CN for a given C-N bond length (in Å).All calculations were performed in the cc-pVDZ basis, energies are listed in a.u., and FC denotes the frozen-core approximation.DMRG(FC) calculations were performed with the DBSS procedure and TRE=10 −5 .