Accurate ab Initio Spin Densities.

We present an approach for the calculation of spin density distributions for molecules that require very large active spaces for a qualitatively correct description of their electronic structure. Our approach is based on the density-matrix renormalization group (DMRG) algorithm to calculate the spin density matrix elements as a basic quantity for the spatially resolved spin density distribution. The spin density matrix elements are directly determined from the second-quantized elementary operators optimized by the DMRG algorithm. As an analytic convergence criterion for the spin density distribution, we employ our recently developed sampling-reconstruction scheme [J. Chem. Phys.2011, 134, 224101] to build an accurate complete-active-space configuration-interaction (CASCI) wave function from the optimized matrix product states. The spin density matrix elements can then also be determined as an expectation value employing the reconstructed wave function expansion. Furthermore, the explicit reconstruction of a CASCI-type wave function provides insight into chemically interesting features of the molecule under study such as the distribution of α and β electrons in terms of Slater determinants, CI coefficients, and natural orbitals. The methodology is applied to an iron nitrosyl complex which we have identified as a challenging system for standard approaches [J. Chem. Theory Comput.2011, 7, 2740].

Introduction

In quantum chemistry, the electronic structure of molecules is described by either ab initio wave function methods or density-functional theory (DFT). For large molecular systems such as transition metal complexes, however, wave function based methods are rarely employed due to the corresponding high computational cost (for counterexamples, see refs (1−5)). Hence, the application of DFT became instrumental in theoretical studies of mechanisms in metal-mediated catalysis.[6−14] Yet, the treatment of open-shell systems[16,18,19] and (near-)degenerate states remains a challenge for DFT.[15] Failures of approximate exchange–correlation density functionals in predicting properties of open-shell systems have been traced to the delocalization error and static correlation error,[17,20] which are rooted in an inappropriate behavior of the energy with respect to fractional charges and fractional spins.[21] In addition to the difficult prediction of ground states from states of different spin,[18,22−30] spin density distributions considerably depend on the approximate exchange–correlation density functional if transition metal complexes containing noninnocent ligands are considered.[31] Qualitatively correct spin density distributions are difficult to obtain within the standard Kohn–Sham formalism that has not been formulated to also produce accurate spin densities.[32]

However, accurate spin densities are desirable for various reasons. (1) In electron paramagnetic resonance spectroscopy (EPR),[33] the spin density is the central quantity on which EPR parameters explicitly depend.[34] Obviously, reliable spin density distributions are important for an accurate calculation of EPR properties, but this remains a difficult task to achieve for theoretical chemistry.[35−39] (2) The question of which approximate exchange–correlation density functional yields sufficiently accurate spin densities remains inconclusive.[31,40] If accurate reference spin density distributions were available, a more detailed analysis of the spin density distribution in terms of spin density difference plots could be used as a qualitative and quantitative benchmark for the validation of approximate exchange–correlation density functionals. (3) According to the Hohenberg–Kohn theorem,[41] the spin density is not needed to calculate the electronic energy or any other expectation value. However, in open-shell systems, it is often introduced as an additional variable which leads to a spin-DFT formalism[42] first introduced by von Barth and Hedin.[43] In spin-DFT, the spin density becomes a fundamental quantity, and reliable reference spin densities could be used to construct proper approximations to the exact exchange–correlation density functional.

For accurate spin densities in cases for which a DFT description fails, ab initio electron correlation methods need to be applied. Pierloot et al. presented complete-active-space self-consistent-field (CASSCF) studies for large transition metal complexes which provided deeper insights into the quality of DFT spin density distributions.[2,44] The large molecular size of these systems requires large active orbital spaces, but the standard CASSCF approach restricts their dimension, which represents the most crucial approximation in such calculations.[40] It is therefore important to understand whether the spin density is converged with respect to the dimension of the active orbital spaces used so far. This is a task that is difficult to study within a standard CAS-type approach.

In general, up to about 18 electrons correlated in 18 spatial orbitals are computationally feasible for standard CASSCF. These limitations may restrict the accurate description of electronic structures which could be approved only by enlarging the dimension of the active orbital space. Reliable reference spin density distributions for complicated open-shell structures as found, for instance, in iron complexes with noninnocent ligands require capabilities beyond those of standard correlation methods.

A different approach for the calculation of correlated ab initio spin densities for large molecules was recently presented by Kossmann and Neese[45] who discussed the performance of orbital-optimized Mo̷ller–Plesset perturbation theory in calculating hyperfine coupling constants for atoms and small molecules. In this approach, isotropic hyperfine constants of coupled-cluster singles-doubles quality could be obtained, which could be further improved by applying spin-component scaling.

Here, we pursue a different route for the calculation of ab initio spin densities by applying the density-matrix renormalization group (DMRG) algorithm. With the DMRG algorithm, introduced by White[46,47] in 1992, much larger active orbital spaces can be considered beyond the limit of, say, 18 electrons correlated in 18 molecular orbitals. It was shown that DMRG is capable of providing accurate wave functions and energies, even for complicated electronic structures (see refs (48−51) for reviews). Moreover, we first showed that the DMRG algorithm yields reliable relative electronic energies between different spin states or isomers of transition metal complexes and clusters for which DMRG was not meant to work and which are a very challenging task for any other multireference quantum chemical method[52] (see also ref (53) for latest results and further references). We shall demonstrate in this work that also accurate DMRG spin density distributions can be determined for very large active orbital spaces.

Recently, we presented a convergence analysis of the spin density distribution for a small iron nitrosyl model complex [Fe(NO)]2+ in a field of point charges, which demonstrated that medium-sized active orbital spaces are sufficient for quantitatively correct spin densities.[40] However, a quantitative analysis that can explore truly large active spaces is still lacking for this complex, which shall therefore be the target system in this work. In such cases, DMRG spin densities can be considered as reliable references which can serve as benchmark results for approximate exchange–correlation density functionals.

This work is organized as follows. In section 2, we discuss the spin density matrix and its spatially resolved counterpart, the spin density distribution, employing the formalism of second quantization. Then, we continue with the introduction of DMRG spin densities. In section 2.2, we present our approach of approximating the DMRG spin density distribution via one from a complete-active-space configuration-interaction(CASCI)-type wave function which allows us to compare DMRG spin densities from calculations with different DMRG parameter sets. In order to validate our approach, we study the spin densities of a medium-sized active orbital space in section 3. This is then extended by considering up to 29 active orbitals in section 4. Finally, a summary and concluding remarks are given in section 5.

Spatially Resolved, Nonrelativistic Spin Densities

Since DMRG is based on the second quantized formalism, we briefly discuss how the spin density in spatial coordinates can be written in second quantization. In first quantization, the operator for the spin density readswhere ŝ is the z component of the one-electron spin operator, r is the spatial coordinate of electron i, and N is the total number of electrons in the system. Applying an orbital basis, the corresponding operator expression in second quantization is given bywhere p,q run over the total orbital basis {ϕ} with ϕ(r) representing the spatial part of a spin orbital. The operators a† and a are the creation and annihilation operators, respectively, for an electron of spin σ in orbital i. In eq 2, the spin density operator δ̂spin(r) is defined in terms of the spin tensor excitation operatorsin the orbital basis (see ref (54) for details). The spatially resolved spin density ρspin(r) is calculated as the expectation value of δ̂spin(r):where |Ψ⟩ represents some normalized reference state|n⟩ = |n1n2...n⟩ is an occupation number vector with elements n ∈ {0,1}. {n} represents the set of all occupation number vectors constructed from k one-particle states. The expectation value on the right-hand side of eq 4 is a spin density matrix element T(:

DMRG Spin Densities from Second-Quantized Elementary Operators

If the reference state |Ψ⟩ is a DMRG wave function in eq 6, the corresponding DMRG spin density matrix elements T(M[DMRG]) are obtained. The matrix representations of the creation and annihilation operators are available in every step of the DMRG algorithm, and each spin density matrix element can thus be easily determined.

The operator a†a in its matrix representation is calculated as a tensor product for which we have to distinguish two different cases. The molecular orbitals p and q are defined either (i) on the same or (ii) on different subsystems of the DMRG partitioning of the active orbital space into the active (sub)system, its environment (the complementary subsystem), and one or two explicitly treated orbitals in between. While the former case is straightforward to handle, for an operator expression in the latter case, however, we need to build operators for the superblock where all three subsystems, i.e., the active subsystem, the exactly represented sites, and the environment, are combined as tensor products.

To illustrate this concept, let us consider two operators a1 and a2 defined on three different subspaces , , and . Then, the combined subspace is defined as = ⊗ ⊗ , where as well as , , and are all subspaces of the N-particle Fock space . For instance, the operator expressions for the combined subspace are given bywhere is the anticommutation matrix of the corresponding subspace . For the product of two operators, we obtainwhere we have used the mixed-product property for the right-hand side of the above equation, which mixes the ordinary matrix product with the tensor product. All remaining operator products can be derived in a similar way. After the spin density matrix is determined, the spatially resolved spin density distribution can be calculated from eq 4. If the wave function is real, the spin density matrix is symmetric and the calculation can be speed up by calculating the upper triangular part of the spin density matrix only.

Spin Density from a Reduced Dimensional CASCI-Type Wave Function

Since CI vectors are in general sparse[55−57]—if contributions below a predefined threshold are neglected—CASCI-type wave functions can be efficiently and accurately projected onto a smaller set of Slater determinants, which only represent the most important contributions to the wave function expansion. We recently reported the sampling-reconstruction algorithm for CASCI-type wave functions defined in a complete active orbital space from a previously optimized DMRG wave function (SRCAS algorithm).[57] An approximate CASCI-type expansion |Ψ̃⟩ for any wave function |Ψ⟩ consisting of k one-particle states can thus be written aswhere the sum runs over all occupation number vectors ñ living in the sampled subspace of the total many-particle Hilbert space. Using eq 6, we can calculate the spin density matrix by substituting the reference state |Ψ⟩ with the approximate state |Ψ̃⟩:

Since the occupation number vectors are orthonormal to each other, the expectation value on the right-hand side of eq 11 can be easily evaluated, and we obtainwhere represents the occupation number vector where orbital p lacks one electron with σ spin. Furthermore, we introduced a phase factor εσ to account for the annihilation operations of a† acting on the bra-state and a acting on the ket-state.

On the basis of this approximate expression for the spin density matrix, we can determine spin density distributions for subspaces of the many-particle Hilbert space of different dimensions and study the sensitivity of the spin density distribution to the number of active-system states in DMRG calculations.

Measures for Spin Density Comparisons

For various reasons, we need suitable measures to assess the similarity of different spin densities. For instance, such a measure would be required to assess the accuracy of a given spin density compared to a reference spin density.

Monitoring the evolution of the spin density for an increasing number of active-system states[58] can illustrate the convergence behavior of the spin density distribution with respect to the number of active-system states m. Isosurface plots of the difference in spin density distributions for calculations with different m values can only serve as a qualitative convergence measure. As quantitative measures, however, we introduce two distances which quantify how far two spin densities are apart from each other. Both distance measures are defined with the absolute error in the spin density difference distribution. The accumulated absolute error Δabs is given byand the root-square error Δrs readswhere ρspin(r) refers to the spin density distribution corresponding to some calculation i, e.g., to a CASSCF or DMRG spin density when different chemical methods are compared, or to some parameter sets if different spin densities are determined with the same method. If two spin densities ρspin(r) and ρspin(r) are similar, both Δabs and Δrs approach zero. For accurate ab initio spin densities, we shall require both error measures to be smaller than 0.005 (Δabs) or 0.001 (Δrs), respectively (in view of the results discussed in section 3.1).

A different similarity measure can be applied by employing directly the knowledge of the reconstructed CASCI-type wave function expansion. This procedure relies on the closeness measure of two quantum states, namely the quantum fidelity.[59,60] The importance and potential application of the quantum fidelity within the DMRG framework was first discussed by some of us[61] in the context of quantum error correction and was also utilized in our SRCAS approach.[57] Two CASCI-type wave function expansions reconstructed for different numbers of DMRG active-system states, m1 and m2, can be explicitly compared by calculating their quantum fidelityas an overlap measure.

A Noninnocent Model System

In a previous study, we reported DFT and CASSCF spin density distributions in iron nitrosyl complexes as well as for the [Fe(NO)]2+ molecule embedded in a square-planar field of point-charges to emulate the one-electron states of the full complexes.[40] Since DFT spin densities of iron nitrosyl complexes remain ambiguous, we choose the small [Fe(NO)]2+ molecule in its doublet state for our analysis here. The point charges facilitate a dynamic change of the character of the electronic wave function by shortening the distances dpc of the point charges to the metal center. Depending on this distance dpc, both single-reference and multireference situations can be created for [Fe(NO)]2+. When the four point charges are located at a distance of 1.131 Å from the iron atom, the electronic structure of the [Fe(NO)]2+ molecule represents a single-reference problem, while for dpc = 0.598 Å, a multireference case is generated.

The [Fe(NO)]2+ structure features a Fe–N bond length of 1.707 Å and a N–O bond distance of 1.177 Å with a Fe–N–O angle of 146°. The four negative point charges of −0.5e each are located as depicted in Figure 1e. Due to the small size of the [Fe(NO)]2+ molecule, we can efficiently study the dependence of the spin density distribution on different DMRG parameter sets such as the number of DMRG active-system states m. Thereby, we are able to define appropriate convergence measures for the spin density in order to reach a predefined accuracy.

Figure 1

(a and c) Spin density difference plots for DMRG(7,7)[m] spin densities calculated for a different number of DMRG active-system states m with respect to the CAS(7,7)SCF reference spin densities shown in (b) and (d), respectively, for [Fe(NO)]2+ in a quadratic-planar point-charge field. Two different distances are considered, namely, dpc = 1.131 Å (a and b) and dpc = 0.598 Å (c and d). For both distances, the spin density could be perfectly reproduced in a DMRG(7,7)[64] calculation and is therefore not depicted here. An isosurface value of 0.0003 in (a) and (b) and 0.003 in (c) and (d), respectively, is chosen. (e) Structure of [Fe(NO)]2+ with the four point charges of −0.5e. In this and in the following figures, the blue (yellow) color denotes an excess of α-electron density, while yellow (blue) represents an excess of β-electron density in the spin density (difference) plots.

The Single-Reference Case

As already discussed in great detail in ref (40), the minimal active orbital space for [Fe(NO)]2+ with dpc = 1.131 Å comprises seven electrons correlated in seven orbitals for qualitatively reliable spin density distributions. It consists of Fe 3d (d, d, d, d and d) and both NO π* orbitals. As an orbital basis in our DMRG calculations, the natural orbitals from a CAS(7,7)SCF calculation performed with the Molpro program package[62] using Dunning’s cc-pVTZ basis set for all atoms[63,64] were taken. The one-electron and two-electron integrals in the natural orbital basis were also calculated with the Molpro program package.[62] All DMRG calculations reported in this section were carried out with the Zurich DMRG program.[65] Random noise was added to the density matrix in order to force the mixing of configurations that would have not been captured otherwise if the number of active-system states m was too small.[66,67]

We performed DMRG calculations for different numbers of DMRG active-system states m abbreviated as DMRG(x,y)[m], where x corresponds to the number of active electrons and y is the number of active orbitals for m renormalized active-system states. Starting with m = 16, m is further increased to 32 and 48 until the CAS(7,7)SCF reference energy is reproduced for m = 64 active-system states (see Table 1). Note that the number of active-system states needed to reproduce the CASSCF result is very small in this case. This can be explained employing concepts of quantum information theory illustrated in section 4. The DMRG calculations reported in this section do not employ these concepts to enforce better convergence. This decision is deliberately made in order to produce nonconverged low-m results to compare with the m = 64 calculation. We should note that this artifact could be cured by the dynamical block state selection (DBSS) procedure,[68,69] while the strong dependence on small m values and the convergence to local minima can be overcome by applying the configuration interaction-based dynamically extended active space (CI-DEAS) procedure.[61]

Table 1

Ground State Energy for [FeNO]2+ Surrounded by Four Point Charges at Two Different Distance Sets dpc in Hartree Atomic Units for CAS(7,7)SCF and DMRG(7,7)[m] Calculations for Different Numbers of DMRG Active-System States m

dpc = 1.131 Å		dpc = 0.598 Å
method	E/Hartree	method	E/Hartree
HF	–1392.844043	HF	–1396.821220
CAS(7,7)SCF	–1392.887247	CAS(7,7)SCF	–1396.858313
DMRG(7,7)[16]	–1392.881067	DMRG(7,7)[16]	–1396.762709
DMRG(7,7)[32]	–1392.885462	DMRG(7,7)[32]	–1396.818651
DMRG(7,7)[48]	–1392.886893	DMRG(7,7)[48]	–1396.840018
DMRG(7,7)[64]	–1392.887247	DMRG(7,7)[64]	–1396.858313

The spin density distributions for our four DMRG calculations (m = 16, 32, 48, 64) are shown in Figure 1a and were determined as discussed in section 2.1. To emphasize the dependence on m, the corresponding spin density difference plots with respect to the CAS(7,7)SCF reference spin density distribution (shown in Figure 1b) are displayed. Note that all isosurface plots are shown for the same isosurface value of 0.0003, where a blue surface corresponds to an excess of α-electron density, while a yellow surface corresponds to an excess of β-electron density for all spin density isosurface plots shown. All DMRG calculations yield qualitatively similar spin density distributions; only minor quantitative differences can be observed. The CAS(7,7)SCF reference spin density can be perfectly well reproduced for m = 64 DMRG active-system states and is, hence, not shown in Figure 1a.

To calculate approximate spin density distributions from reconstructed CASCI-type wave functions, we first have to sample the most important configurations of the N-particle Hilbert space. For this purpose, we applied our SRCAS method.[57] Due to the small size of the active space, the N-particle Hilbert space is spanned by only 1225 Slater determinants, and all corresponding CI coefficients can be determined directly from the CASSCF reference calculation. In general, similar CI coefficients are obtained for all DMRG calculations and the CASSCF reference; i.e., similar wave functions are converged, resulting in small differences in the spin density distributions. The distribution of the CI coefficients is depicted in Figure 1 of the Supporting Information.

Spin density distributions determined for different sampled subspaces of the N-particle Hilbert space are in good agreement with the corresponding DMRG spin density. Note that the sampled subspaces are defined by the threshold value of the completeness measure (COM) introduced in ref (57) with COM = (1 – ΣC2), where I runs over all sampled configurations with CI coefficients C. In general, threshold values of 0.01 to 0.001 turned out to be sufficient for obtaining quantitatively reliable spin densities in this single-reference case. The corresponding isosurface plots and excitation histograms with respect to the COM are summarized in the Supporting Information.

The spin density difference plots in Figure 1a illustrate the convergence of the spin density distribution with respect to the number of DMRG active-system states m. The absolute error Δabs and the root-square error Δrs of the spin density difference distributions provide a quantitive measure for the accuracy (see Table 2). The differences in the spin densities calculated for m = 48 DMRG active-system states is small compared to the CAS(7,7)SCF reference. For 48 active-system states upward, both Δabs and Δrs are below their threshold values, given in section 2.3. The set of quantum fidelity measures F for our four DMRG calculations with m ∈ {16, 32, 48, 64} is {0.980000, 0.994395, 0.999012}. Increasing m from 48 to 64 DMRG active-system states corresponds to F48,64 = 0.999012, which illustrates the similarity of both DMRG wave functions and results in reliable spin density distributions for m ≥ 48.

Table 2

The Absolute Error Δabs and the Root-Square Error Δrs of the DMRG(7,7)[m] Spin Densities with Respect to the CAS(7,7)SCF Reference for [FeNO]2+ Surrounded by Four Point Charges at Two Different Distance Sets dpc Employing Different Numbers of DMRG Active-System States m

	dpc = 1.131 Å		dpc = 0.598 Å
method	Δabs	Δrs	Δabs	Δrs
DMRG(7,7)[16]	0.007678	0.002147	0.213543	0.052168
DMRG(7,7)[32]	0.004392	0.001285	0.221198	0.052144
DMRG(7,7)[48]	0.001397	0.000412	0.081631	0.020418
DMRG(7,7)[64]	1.34 × 10–5	5.40 × 10–6	9.69 × 10–6	3.66 × 10–6

The Multireference Case

A multireference character of the [Fe(NO)]2+ molecule can be induced by decreasing the distances of the point charges to the iron atom. In the squeezed model complex, the point charges are placed at a distance of dpc = 0.598 Å from the iron center in the same configuration as before. Similar to the single-reference problem, the minimum active orbital space considered here comprises seven electrons correlated in seven orbitals. Yet, it consists of four Fe d, d, d, and d (d is excluded due to the compressed point charge environment), two NO π*, and one NO σ orbital which interacts with the Fe d orbital. Again, the natural orbitals from a CAS(7,7)SCF calculation were taken as orbital bases in our DMRG calculations and determined with the Molpro program package[62] using Dunning’s cc-pVTZ basis set for all atoms[63,64] The calculation of the one-electron and two-electron integrals in this natural orbital basis was also performed with the Molpro program package.[62] All DMRG calculations were carried out with the Zurich DMRG program.[65] As before, we performed DMRG calculations for four different numbers of DMRG active-system states m. Starting with m = 16, m is further increased to 32 and 48 until the CAS(7,7)SCF reference energy is obtained for m = 64 active-system states (see Table 1). The small-m calculations are designed not to reproduce the CAS(7,7)SCF reference for this analysis. Note, however, that a small number of active-system states was sufficient to reproduce the CASSCF result as observed in the single-reference problem.

In Figure 1d, the CAS(7,7)SCF spin density distribution is shown, which is taken as the reference distribution, while Figure 1c illustrates the spatially resolved differences in the DMRG(7,7)[m] and CAS(7,7)SCF spin density distributions. Note that the same isosurface value of 0.003 was chosen for all spin densities shown. For small m values, qualitatively different spin density distributions are obtained. The β-electron density around the nitrosyl ligand is underestimated and a dumbbell-shaped β-electron density is obtained in contrast to the cylindric shape of the reference β-electron density. The α-electron density around the Fe atom is underestimated. Increasing m to 48 results in a cylindric β-electron density around the NO ligand which differs only little from the reference spin density. The spin density can be exactly reproduced for m = 64 active-system state for which also the CAS(7,7)SCF reference energy is obtained. The convergence properties of the DMRG(7,7)[m] spin density with respect to m can be quantified by the Δabs and Δrs values where significantly large values (>0.005 and >0.001, respectively) are obtained for spin density distributions determined in small-m calculations (see Table 2).

In Figure 2, the distribution of CI coefficients for the DMRG and CASSCF wave functions is shown. Since only the position of the point charges has been modified, the N-particle Hilbert space remains spanned by 1225 Slater determinants, and all corresponding CI coefficients can be determined directly from the CASSCF reference calculation as in the single-reference case. Similar CI coefficients are obtained for the DMRG(7,7)[64] calculation and the CAS(7,7)SCF reference; i.e., similar wave functions are converged. However, significantly different CI coefficients are obtained—as expected—for smaller m values. In particular, the deviations are most significant for configurations corresponding to the largest (absolute) CI weights. Additional information on the distribution of CI coefficients can be found in the Supporting Information.

Figure 2

Distribution of the absolute value of the CI coefficients corresponding to the Slater determinants in the DMRG(7,7)[m] calculations with different renormalized active-system states m and in the CAS(7,7)SCF reference calculation for [Fe(NO)]2+ surrounded by four point charges at a distance of dpc = 0.598 Å from the iron atom. All Slater determinants are ordered according to the CI weights of the CAS(7,7)SCF calculation.

Similarly to the single-reference problem discussed above, reliable spin densities obtained from reduced dimensional CASCI-type wave function expansions can be determined for a COM ≥ 0.001 independent of m. A complete collection of spin density distributions for different CASCI-type wave function expansions and DMRG parameter sets can be found in the Supporting Information. Figure 3 shows the ratio of Slater determinants with respect to the complete N-particle Hilbert space which have been picked up in the sampling procedure and sorted by their corresponding CI weights for the DMRG(7,7)[32] and DMRG(7,7)[64] calculation. For COM ≥ 0.001, the reconstructed CASCI-type wave function contains the major part of the important Slater determinants, while for a further decreased threshold value of 10–5 almost all significant Slater determinants have been picked up. Note that the sampling procedure was restricted to accept only configurations with (absolute) CI coefficients larger than the threshold value for COM. Although all possible excitations are included in the CASCI-type wave function in the limit of COM → 0 (see also Figure 4 of the Supporting Information), the pattern of the CI coefficients of the DMRG(7,7)[64] calculation is different from the CI pattern of the CAS(7,7)SCF reference. While large CI coefficients (|C| > 0.0001) are reproduced within sufficient accuracy, smaller CI weights are underestimated. The maximum of the curve is shifted toward smaller CI weights <10–7. Hence, the DMRG algorithm disregards an exact weighting of unimportant configurations with small CI coefficients, which is a feature of matrix product and tensor network states where large CI coefficients should be reproduced and unimportant configurations are neglected[70,71] A complete collection of excitation histograms for different CASCI-type wave functions can be found in the Supporting Information. To quantify the differences in the underlying wave functions for our four DMRG calculations employing m ∈ {16, 32, 48, 64} active-system states, we calculated the quantum fidelity F, which forms in this case a set of overlap measures of {0.831887, 0.897445, 0.955669}.

Figure 3

CI histogram of the absolute values of the CI coefficients for the Slater determinants for reconstructed CASCI-type wave function expansions from the DMRG(7,7)[m] calculations with different renormalized active-system states m for the [Fe(NO)]2+ molecule surrounded by four point charges at a distance of dpc = 0.598 Å from the iron atom. The CAS(7,7)SCF reference calculation is also shown for comparison. thr corresponds to the threshold value of COM in the sampling-reconstruction procedure and denotes the accuracy of the reconstructed CASCI-type wave function. All Slater determinants with CI coefficients in an interval as indicated on the abscissa are grouped together.

Figure 4

Mutual information and single orbital entropies s(1) for the DMRG(13,y)[64] calculations determined for different numbers of active orbitals in the [Fe(NO)]2+ molecule surrounded by four point charges at a distance of dpc = 1.131 Å from the iron center.

We conclude that reliable spin density distributions can be calculated either from converged DMRG ground state wave functions or from the reconstructed CASCI-type wave function expansions. In particular, a fully converged DMRG wave function is not mandatory to obtain qualitatively correct spin density distributions if the CI weights of the most important configurations are well reproduced for a given m value. This holds for both the single-reference and the multireference case. A representative set of Slater determinants, i.e., the most important ones (|C| > 0.001), is sufficient for a qualitatively correct spin density distribution.

Spin Density Distributions for Large Active Spaces

While we have studied the convergence features of DMRG calculations for small active spaces, for which we could obtain an exact CASSCF reference result, we shall now proceed to explore territory with DMRG that is not accessible to the CASSCF approach. In our recent analysis of CASSCF spin densities for the [Fe(NO)]2+ molecule,[40] the spin density distribution was qualitatively converged with respect to the dimension of the active orbital space. For quantitatively accurate spin densities, we need to increase the dimension of the active orbital space so that important iron and ligand orbitals which are missing in the standard CASSCF calculations, e.g., the Fe d double-shell orbital, could also be included in the active orbital space. Here, we extend the convergence series presented in ref (40) by considering active orbital spaces containing up to 29 active orbitals. Starting with an active orbital space comprising 13 active electrons correlated in 20 active orbitals, the number of active orbitals is further increased to 24 and 29, respectively. The two largest active orbital spaces do also contain the fifth d-double-shell orbital which could not be included in all CASSCF calculations presented in ref (40). The [Fe(NO)]2+ molecular structure features the same bond distances and angles as presented in section 3. The four point charges of −0.5e are located at a distance of 1.131 Å from the metal center in order to properly model the square-planar ligand field of the full-fledged complexes in a doublet spin state.

For all DMRG calculations, the natural orbitals from a CAS(11,14)SCF calculation are employed as orbital bases.[72−74] Similarly, the CASSCF calculation as well as the calculations of the one-electron and two-electron integrals in the natural orbital basis were performed with the Molpro program package[62] using Dunning’s cc-pVTZ basis set for all atoms,[63,64] while the DMRG calculations are performed with the Budapest DMRG program.[75] In addition, the DMRG orbital orderings were optimized for all three active orbital spaces, and the CI-DEAS starting guess was performed. Figure 4 displays the corresponding single orbital entropies given byand mutual information determined bywhere i = 1...k is the orbital index and runs over all k one-particle states and ωα, is the α eigenvalue of the reduced density matrix of orbital i,[61] while s(2) is the two-orbital entropy between a pair (i,j) of sites introduced by Rissler et al. to the quantum chemical DMRG algorithm.[76] Note that the mutual information and single orbital entropies are confined to the first 10 natural orbitals for all considered dimensions of the active orbital space. These natural orbitals are highly entangled and represent the most important orbitals comprised in the active orbital space. Therefore, accurate DMRG spin densities can be obtained already for a reasonably small number of active-system states. Similar entropy profiles can be obtained for smaller dimensions of the active orbital space.

The number of DMRG active-system states m was set to 128, 256, 512, 1024, and 2048, respectively. The ground state energies for all DMRG calculations are summarized in Table 3. Considering the DMRG(13,20)[m] calculations, an energy convergence of 0.135 mH (0.4 kJ/mol) is reached with respect to m. For the largest active orbital space, the DMRG(13,29)[1024] energy is converged to 1.195 mH (3.1 kJ/mol) when compared to the DMRG(13,29)[2048] reference.

Table 3

Ground State Energy for [Fe(NO)]2+ Surrounded by Four Point Charges at a Distance of dpc = 1.131 Å from the Iron Center in Hartree Atomic Units for Our DMRG(x,y)[m] Calculations Employing Different Numbers of DMRG Active-System States ma

E/Hartree
method	DMRG(13,20)	DMRG(13,24)	DMRG(13,29)
m = 128	–1393.014662	–1392.991085	–1393.014010
m = 256	–1393.018626	–1393.019309	–1393.024883
m = 512	–1393.020065	–1393.021876	–1393.030374
m = 1024	–1393.020511	–1393.022946	–1393.033001
m = 2048	–1393.020646	–1393.023294	–1393.034196

The CAS(11,14)SCF energy is −1393.013 396 Hartree.

Convergence of DMRG Spin Densities

The dependence of the spin density distribution on the number of DMRG active-system states m is shown in Figure 5 where the differences in spin density distribution are plotted for DMRG(13,y)[m] calculations with respect to the converged DMRG(13,29)[2048] reference calculation. For increasing m values, the differences in the spin density distribution decrease (see each row in Figure 5 from the left to the right). Similarly, we observe that the spin density gradually converges with respect to the dimension of the active orbital space (see last column from the top to the bottom of Figure 5). In particular, changes in the spin density are negligible when m is increased from 1024 to 2048, and hence, reliable spin density distributions can be obtained even if the total energy is not yet converged with respect to m (the difference is 1.195 mH, see above).

Figure 5

DMRG(13,y)[m] and CAS(x,y)SCF spin density difference plots with respect to the DMRG(13,29)[2048] spin density distribution (d) for [FeNO]2+ surrounded by four point charges at a distance of dpc = 1.131 Å from the iron center. All spin densities are displayed for an isosurface value of 0.001. (a) DMRG(13,20)[m]–DMRG(13,29)[2048] spin density difference plots. (b) DMRG(13,24)[m]–DMRG(13,29)[2048] spin density difference plots. (c) DMRG(13,29)[m]–DMRG(13,29)[2048] spin density difference plots. (d) The DMRG(13,29)[2048] reference spin density distribution.

Furthermore, the Δabs and Δrs values quantify the convergence series of the determined DMRG spin density distributions. In Table 4, both error quantities are listed for each DMRG(x,y)[m] spin density with respect to the DMRG(13,29)[2048] reference spin density. In general, the absolute error Δabs and the root-square error Δrs decrease for increasing m, keeping the dimension of the active orbital space fixed. Note that larger active orbital spaces require a larger m value to obtain the same accuracy as achieved in smaller active space calculations. This is not immediately evident from the error data presented in Table 4 since different dimensions of the active orbital space are compared, which result in nonzero error values, while error values determined for different parameter sets, but the same dimension of the active orbital space, could vanish. The large error values for the DMRG(13,24)[128] calculation indicate that important states were not picked up by the DMRG algorithm, resulting in the large differences in the spin density distribution displayed in Figure 5. Furthermore, since both error values determined for the DMRG(13,29)[1024] calculation are below the threshold values, no considerable improvement in the accuracy of the spin density distribution can be expected when m is further increased to more than 2048 active-system states.

Table 4

The Absolute Error Δabs and the Root-Square Error Δrs of the DMRG(13,y)[m] Spin Densities with Respect to the Converged DMRG(13,29)[2048] Reference Spin Density for a Different Number of Normalized Active-System States m for [FeNO]2+ Surrounded by Four Point Charges at a Distance of dpc = 1.131 Å from the Iron Centera

method	Δabs	Δrs
DMRG(13,20)[128]	0.030642	0.008660
DMRG(13,20)[256]	0.020088	0.004930
DMRG(13,20)[512]	0.016415	0.003564
DMRG(13,20)[1024]	0.015028	0.003162
DMRG(13,20)[2048]	0.014528	0.003028
DMRG(13,24)[128]	0.590022	0.235922
DMRG(13,24)[256]	0.020993	0.003245
DMRG(13,24)[512]	0.014045	0.003633
DMRG(13,24)[1024]	0.011622	0.002668
DMRG(13,24)[2048]	0.010731	0.002361
DMRG(13,29)[128]	0.032171	0.010677
DMRG(13,29)[256]	0.026005	0.006790
DMRG(13,29)[512]	0.010826	0.003406
DMRG(13,29)[1024]	0.003381	0.000975
CAS(11,11)SCF	0.086658	0.024495
CAS(11,12)SCF	0.080249	0.020591
CAS(11,13)SCF	0.046303	0.011402
CAS(11,14)SCF	0.042544	0.010954
CAS(13,13)SCF	0.052239	0.012124
CAS(13,14)SCF	0.073400	0.019850
CAS(13,15)SCF	0.053157	0.011180
CAS(13,16)SCF	0.104928	0.031922

The Δabs and Δrs values of the CAS(x,y)SCF calculations of ref (40) with respect to the DMRG(13,29)[2048] reference spin density are also listed.

In order to demonstrate the convergence of the DMRG(13,29) wave function with respect to m (and thus the convergence of the obtained DMRG(13,29)[2048] reference spin density distribution), the CASCI-type wave function expansions are reconstructed and compared for all m values. In particular, the influence of the missing d-double-shell orbital can be assessed by examining the CI coefficients corresponding to Slater determinants with an occupied d-double-shell orbital. Following the conclusions of a benchmark study for intermediate CAS sizes (see Supporting Information), only the most important configurations (|C| ≥ 0.00001) are necessary to obtain an accurate wave function expansion. As the convergence threshold for the sampling procedure, a value of 0.001 is sufficient. With this threshold, similar CASCI-type wave function expansions are obtained for a quantum fidelity measure close to 0.998. The set of quantum fidelity measures F for our five DMRG calculations with m ∈ {128, 256, 512, 1024, 2048} is {0.991800, 0.995510, 0.996983, 0.997639}. As the number of DMRG active-system states is enlarged, the CI coefficients of the reconstructed wave function expansion converge gradually, which is indicated by the increasing quantum fidelity measure. Note that F is close to the ideal value of 0.998 already for a small number of DMRG active-system states m, and hence, only minor variations in the large CI coefficients occur when m is increased, which explains the slight differences in the spin density distributions displayed in Figure 5c.

To demonstrate that this is indeed the case, the CI coefficients of the most important Slater determinants (|C| > 0.0001) corresponding to the m = 128 and m = 1024 calculations are shown in Figure 6. Slater determinants with large CI weights (|C| > 0.05) are similar for both DMRG parameter sets; only minor deviations can be observed. Note that all of these Slater determinants have been incorporated in the DMRG wave function already for m = 128. Considerable differences in CI weights are present for Slater determinants corresponding to small-valued CI coefficients (|C| < 0.015), while some Slater determinants with |C| < 0.01 have not been incorporated in the DMRG wave function for m = 128. These off-size or missing configurations lead to the different spin density distributions for small m values.

Figure 6

Distribution of the absolute value of the CI coefficients for the DMRG(13,29)[m] calculations with m = 128 and 1024, respectively, for [FeNO]2+ surrounded by four point charges at a distance of dpc = 1.131 Å from the iron center. The CI coefficients reconstructed for both DMRG calculations are always printed for the same Slater determinants. The determinants are ordered according to the CI weight of the DMRG(13,29)[2048] reference calculation.

From the reconstructed CASCI-type wave function, the influence of the d-double-shell orbital as well as of the empty ligand orbitals on the spin density distribution can be analyzed. In the upper part of Table 5, configurations containing an occupied d-double-shell orbital and corresponding to the largest CI coefficients are presented. In the lower part of Table 5, some selected configurations with large CI coefficients carrying excitations to empty ligand orbitals that cannot be included in standard CASSCF calculations are presented and compared for the DMRG(13,29)[128] and DMRG(13,29)[1024] calculations. In general, Slater determinants with an occupied d-double-shell orbital feature small CI weights (|C| ≤ 0.003) and are hence of minor importance, while Slater determinants bearing occupied ligand orbitals feature large CI coefficients. Configurations containing occupied ligand orbitals that are only included in the DMRG(13,29)[m] calculations (marked in bold face in Table 5) possess considerably large CI weights. All other Slater determinants with excitations to different empty ligand orbitals have smaller CI coefficients. Hence, those ligand orbitals pose a significant contribution in obtaining accurate spin density distributions for the small model complex and cannot be neglected from the active orbital space.

Table 5

Some Important Occupation Number Vectors (ONV) with the Corresponding CI Weights from DMRG(13,29)[m] Calculations for [FeNO]2+ Surrounded by Four Point Charges at a Distance of dpc = 1.131 Å from the Iron Centera

	CI weight
Slater determinant	m = 128	m = 1024
b2b222a0a0000000 0000000 a 00000	0.003252	0.003991
bb2222aa00000000 0000000 a 00000	–0.003226	–0.003611
222220ab00000000 0000000 a 00000	–0.002762	–0.003328
ba2222ab00000000 0000000 a 00000	0.002573	0.003022
b2a222a0b0000000 0000000 a 00000	–0.002487	–0.003017
202222ab00000000 0000000 a 00000	0.002405	0.002716
b222a2a0b0000000 0000000 0 0000a	0.010360	0.011558
22b2a2a0a0000000 0000000 0 b0000	0.009849	0.011366
22b2a2a0b0000000 0000000 0 a0000	–0.009532	–0.011457
b2222aab00000000 0000000 0 0000a	–0.009490	–0.010991
a2222baa00000000 0000000 0 0000b	–0.009014	–0.010017
b2b222a0a0000000 0000000 0 0a000	0.008820	0.010327
b2222aab00000000 00a0000 0 00000	–0.004277	–0.005436
22b2a2a0b0000000 a000000 0 00000	–0.004224	–0.006852

Upper part: ONVs containing an occupied d double-shell orbital (marked in bold face). Bottom part: additional selected important configurations with occupied natural orbitals that cannot be included in the active orbital space in CASSCF calculations (marked in bold face), for the same DMRG(13,29)[m] calculations. 2: doubly occupied natural orbital. a: natural orbital occupied by an α electron. b: natural orbital occupied by a β electron. 0: empty natural orbital.

Assessment of CASSCF Spin Densities

The converged DMRG(13,29)[2048] reference spin density can be used to assess the accuracy of CASSCF spin density distributions and benchmark the quality of the (restricted) active orbital spaces in standard CASSCF calculations (see Figure 7a). Note that the same isosurface value has been taken to display the DMRG(13,y)[m]–DMRG(13,29)[2048] and CAS(x,y)SCF–DMRG(13,29)[2048] spin density difference plots. The CASSCF spin density distributions determined for medium-sized active orbital spaces oscillate around the converged DMRG spin density. Depending on which double-d-shell orbital is included in the active orbital space, the β-electron density around the NO ligand is either overestimated or underestimated. This results either in pure spin-polarized cases with β-electron density found only around the nitrosyl ligand for CAS(11,11), CAS(11,14), CAS(13,13), and CAS(13,14) or some additional α-electron density present around the NO ligand associated with a simultaneous decrease in the β-electron density for CAS(11,12), CAS(11,13), CAS(13,15), and CAS(13,16).

Figure 7

(a) CAS(x,y)SCF and (b) DFT spin density difference plots with respect to the DMRG(13,29)[2048] spin density distribution for [FeNO]2+ surrounded by four point charges at a distance of dpc = 1.131 Å from the iron center. All spin densities are displayed for an isosurface value of 0.001.

Similarly, the large Δabs and Δrs values stress the differences in the spin density distributions which are considerably larger than those from the DMRG(13,y)[m]–DMRG(13,29)[2048] difference analysis (Table 4). Furthermore, Table 4 indicates that the CAS(11,11)SCF and CAS(11,12)SCF calculations and the CAS(11,13)SCF and CAS(11,14)SCF calculations, respectively, are of similar accuracy, as they have similar error values, but the spin density difference plots emphasize the qualitatively different spin density distributions. Increasing the dimension of the active orbital space results in even larger deviations from the DMRG reference spin density because the active space is not stable and important orbitals are rotated out of the CAS. Note that all DMRG calculations—except DMRG(13,24)[128]—yield smaller error values and smaller differences in the spin density difference plots.

Although the CASSCF spin densities are quantitatively converged with respect to the active orbital space, significant qualitative—but also non-negligible quantitative—differences to the DMRG(13,29)[2048] reference spin density can be observed. The extension of the active orbital space by including an additional shell of d orbitals only is not sufficient to obtain a qualitatively accurate spin density distribution for the small iron nitrosyl molecule. Our analysis indicates that empty ligand orbitals are essential for calculating reliable reference spin densities. This may have severe implications for the standard CASSCF approach that require further analysis in future work.

Comparison to DFT Spin Densities

A comparison of DFT and CASSCF spin density distributions for medium-sized active orbital spaces for the [FeNO]2+ molecule has already been discussed in our previous work (see ref (40) for more details). For an unambiguous benchmark of approximate exchange–correlation density functionals, the DFT spin densities of ref (40) can be compared to the DMRG reference distribution. The qualitative analysis of the DFT–DMRG(13,29)[2048] spin density difference distributions is shown in Figure 7b. When comparing to the results obtained in ref (40), similar conclusions concerning the performance of approximate exchange–correlation density functionals can be drawn. The best agreement is found for BP86, BLYP, and TPSS, while the remaining approximate exchange–correlation density functionals yield larger deviations and result in too large spin polarization. We should note that BP86, BLYP, and TPSS correctly predict the distribution of the α-electron density around the nitrosyl ligand, although it is overemphasized. In general, nonhybrid functionals yield spin densities which are in closest agreement with the DMRG reference distributions. This observation is supported by both error measures which are smallest for BP86, BLYP, and TPSS (see Table 1 in the Supporting Information).

Conclusions and Outlook

In this work, we have demonstrated how reliable ab initio spin density distributions can be calculated for very large active spaces. Our procedure is based on the DMRG algorithm and on two different approaches to obtain spin density matrix elements: (i) on-the-fly directly from the second-quantized DMRG elementary operators or (ii) from an approximate CASCI-type wave function expansion which is determined by our SRCAS algorithm.[57] The reconstructed CASCI-type wave function can also be used as a means to compare a series of DMRG calculations employing a different number of DMRG active-system states m.

The small noninnocent molecule [FeNO]2+ surrounded by four point charges represents a suitable system to validate our approach. The spin density distributions are highly sensitive to the nature of the converged state. We deliberately converged DMRG wave functions that correspond to local minima in the electronic energy in order to compare with qualitatively wrong wave functions. The possibility of convergence into local minima is shown by examining the (largest) CI coefficients of the SRCAS-reconstructed CASCI-type wave function. Strong deviations with respect to the absolute value of the CI coefficients indicate that the number of DMRG active-system states m is chosen too small, and hence important states have not been incorporated by the DMRG algorithm. Spin densities corresponding to such local minima deviate considerably from the ground state spin density.

The convergence analysis of the spin density distribution for the [FeNO]2+ molecule considered active orbital spaces comprising up to 29 active orbitals. Difference plots of the spin density distribution for different active orbital spaces as well as the absolute error and the root-square error in the spin density difference distribution indicate a quantitatively converged spin density with respect to the dimension of the active orbital space and the number of active-system states m (which was as large as m = 2048). The DMRG reference spin density has been used to validate CASSCF spin densities resulting in significant quantitative and even qualitative differences. Considering an additional shell of d orbitals is not sufficient to obtain reliable spin densities for the small model system, and the active orbital space must be extended by additional unoccupied ligand orbitals. Similar difficulties are likely to be present for larger iron nitrosyl complexes where the point charges are replaced by different ligands, and hence additional ligand and iron orbitals must be included in the active orbital space. The DMRG study of larger {FeNO}7 complexes is now pursued in our laboratory.

A convergence analysis of the spin density in terms of spin density difference plots with respect to the number of DMRG active-system states indicates that reliable reference spin densities can be obtained even if total energies are not converged with respect to m. A similar conclusion was found in our previous work regarding the energy splittings of states of different spin multiplicity.[50,52,53] Comparison of CI weights corresponding to the most important configurations of the reconstructed CASCI-type wave functions for different m values furthermore ensures that reliable spin densities are obtained. The similarities in DMRG wave functions can be quantified by the quantum fidelity measure, which can be used as an additional convergence criterion for spin density distributions in a sequence of DMRG calculations.

Spin densities calculated from approximate CASCI-type wave functions are in good agreement with the DMRG reference spin density. Qualitatively reliable spin densities can be obtained even for large thresholds of COM (0.001) when the most important configurations have been picked up in the wave function expansion. For this threshold, the CASCI-type wave function contains Slater determinants with absolute CI weights larger than 0.00001 which are important for the spin density.

The comparison of DFT spin densities with the DMRG reference distributions allows us to benchmark approximate exchange–correlation density functionals. Although nonhybrid functionals yield spin density distributions closest to the DMRG reference, significant qualitative and quantitative differences to the DMRG reference distributions could be observed for all investigated density functionals. Similar conclusions were drawn in our previous study, where DFT spin densities were assessed against CASSCF spin densities,[40] entailing that none of the investigated exchange–correlation density functionals yields sufficiently accurate spin density distributions for the [FeNO]2+ molecule.