Chemical insights into the electronic structure of Fe(II) porphyrin using FCIQMC, DMRG, and generalized active spaces

Stochastic-CASSCF and DMRG procedures have been utilized to quantify the role of the electron correlation mechanisms that in a Fe-porphyrin model system are responsible for the diﬀerential stabilization of the 3 E g over the 5 A 1 g state. Orbital entanglement diagrams and CI coeﬃcients of the wave function in a localised orbital basis allow for an eﬀective interpretation of the role of charge-transfer conﬁgurations. A preliminary version of the Stochastic Generalized Active Space Self-Consistent Field method has been developed and is here introduced to further assess the π backdonation stabilizing eﬀect. By the new method excitations between metal and ligand orbitals can selectively be removed from the complete CI expansion. It is demonstrated that these excitations are key to the diﬀerential stabilization of the triplet, eﬀectively leading to a quantitative measure of the correlation enhanced π backdonation.


Introduction
Fe porphyrins are central building blocks for various enzymes in biochemistry, playing there an important role in redox reactions, and charge and molecular transport [1,2,3,4].Their enzymatic activity is bound to the presence of nearly degenerate electronic states.For example, Fe(II) porphyrin has a d 6 electronic configuration with high-spin (quintet), intermediate-spin (triplet), and low-spin (singlet) states close in energy.The relative stability of these states depends on ligand field and many-body correlation effects experienced by the metal center, that, in turn, depend on chemical functionalization and geometry of the conjugated macrocycle.
In a recent work [5], the Stochastic-CASSCF method [6] has been applied to a simplified Fe(II) porphyrin model complex.Valence 3d, nitrogen σ, and the entire conjugated π system electrons have been explicitly correlated in the Stochastic-CASSCF (32,34) procedure, while variationally optimizing the molecular orbitals.A correlation enhanced σ donation/π backdonation mechanism has been proposed to explain the relative stability of the intermediate-spin state (triplet) over the high-spin (quintet).According to this mechanism, the electron repulsion at the metal center is differentially reduced for the triplet spin state by electron delocalization induced by charge-transfer (CT) excitations between the 3d orbitals at the metal center and the π system of the macrocycle.This mechanism is not captured in smaller active space calculations and it is not accounted for by second order perturbation theory procedures (such as CASPT2) on top of small active space reference wave functions.The correlation enhanced charge transfer mechanism suggests that large active spaces, such as the CAS (32,34) used in reference 5, are crucial for an accurate description of electron delocalization in metal porphyrins.
In two subsequent works [7,8] the role of semi-core correlation (3s3p shell) has been investigated.Phung et al. [7] tackled semi-core correlation effects via a composite CASPT2/CC method, while Li Manni and Alavi [8] further enlarged the valence (32,34) active space, including (3s, 3p) orbitals and their electrons, CAS (40,38).In both cases, inclusion of semi-core correlation had the effect to further stabilize the intermediate spin state over the high spin state.
In this report we analyze in greater detail the electronic structure of the quintet 5 A 1g and triplet 3 E g states to quantitatively assess the electron correlation mechanisms that lead to the stabilization of the triplet over the quintet spin state.
In this work, the optimized CASSCF (32,34) wave functions from Ref. 5 have been represented in the one-electron basis of localized orbitals, allowing to classify the leading electronic configurations of the CAS wave functions into on-site (local) excitations and CT excitations, and connect quantum chemical findings to the chemical concepts of σ donation and π backdonation.Due to the size of the active space and of the corresponding CI expansion the Full Configuration Quantum Monte Carlo (FCIQMC) and the Density Matrix Renormalization Group (DMRG) algorithms have been used as CI eigensolvers, here referred to as Stochastic-CASCI [6] and DMRG-CASCI.
Orbital entanglement measures [9,10,11,12,13,14] based on DMRG for quantum chemistry [15,16,17,18,19,20,21,22,23,24,25] provide a clear description of correlation effects that differentially stabilize the 3 E g over the 5 A 1g spin state.Furthermore, in order to directly relate the role of the CT configurations to the triplet-quintet spin-gap, a strategy has been designed to selectively allow or forbid ligand-to-metal and metal-to-ligand CT excitations within the CAS (32,34) configuration interaction expansion, while letting the rest of the wave function relax within the FCIQMC method [26,27] and orbitals relax self-consistently.This protocol has been inspired by the generalized active space (GAS) procedure with disconnected -GAS subspaces [28], and is to be considered a preliminary version of the Stochastic-GASSCF method.
These restrictions in the active space wave function directly assess specific electron-correlation mechanisms, for example π backdonation, while retaining effects related to mean-field interactions.
The DMRG and FCIQMC methods produce very similar and accurate energy estimates, while providing complementary information towards the understanding of the electronic structure of the investigated model system.On one hand, orbital entanglement measures emerge at no extra computational cost from the DMRG procedure, but would require the additional and expensive sampling of higher order reduced density matrices within the FCIQMC framework.On the other hand, FCIQMC gives direct access to the dominant electronic configurations that characterize the CI wave function of the spin states studied in this work.This would require an extra computationally demanding step in DMRG.[29,30,31] In addition, FCIQMC in the stochastic-GASSCF method, allows to enforce restrictions on the type of excitations with ease, which is not straightforward in DMRG.As mentioned above, inclusion of excitation restrictions allows us to quantify the effect of the restricted excitations on the wave function energy and thus their role in stabilizing the triplet spin state.
The remainder of the article is organized as follows.In Section 2 we discuss the important elements of the FCIQMC and DMRG methods, the applied basis sets, and the novel algorithm for Generalized Active Space wave functions in the context of FCIQMC.Results are presented in Section 3 and discussed in Section 4. A summary of the paper and conclusions can be found in Section 5.

Computational Details
Geometry.All calculations were performed on the same Fe porphyrin model complex as in our earlier work.[5] The pyrrole β-carbon atoms were removed and remaining bonds were capped with H-atoms.The Fe−N bond length was kept at 1.989 Å which is closer to the triplet's equilibrium geometry.The molecule was placed into the xy plane with ligand N-atoms on the x and y axes.The point group of the planar molecule is D 4h , however, calculations in the natural orbital basis were carried out using the D 2h Abelian point group, the largest our programs can handle.In the localised basis all calculations have been performed using the C 1 point group symmetry.Basis Set.Basis set and integral specifications have been chosen as in our earlier work [5].Generally contracted atomic natural orbitals (ANO-RCC) of split-valence triple-ζ plus polarization quality, were utilized (ANO-RCC-VTZP) [32,33], resulting in a total of 707 basis functions.Scalar relativistic effects were included with the Douglas-Kroll-Hess integral correction.Evaluation of electron repulsion integrals has been accelerated by using the resolutionof-identity Cholesky decomposition with a threshold of 1 × 10 −4 E h .[34,35,36,37,38] Active space.The active space consists of nine doubly-occupied π-orbitals, seven empty π * -orbitals, five metal centered 3d-orbitals, four doubly-occupied σ N -orbitals, four empty orbitals of the (4s4p) shell, and five empty d -orbitals (double-shell orbitals).With the d 6 -configuration of the Fe(II) metal center, this implies a (32,34) active space.The CASSCF (32,34) optimized active natural orbitals for the quintet and triplet spin states are depicted in Figure 1 and Figure 2, respectivey.
Localisation procedures.The natural active orbitals from Reference 5 have been localised using the method introduced by Pipek and Mezey [39].Furthermore, Procrustes orthogonal transformations [40,41] of the metal centered active orbitals were performed to make them as similar to the atomic orbital basis as possible.The resulting orbitals for the quintet state are displayed in Figure 3.The triplet localized orbitals are very similar and are reported in the supplementary material.The localised orbital basis neither diagonalize the one body density matrix, so it is not possible to assign occupation numbers, nor do they diagonalize the active space Fock matrix, so it is also not possible to assign orbital energies.Since the CASSCF energy expression is invariant under unitary transformations among the active orbitals, and the orbitals utilized were already the optimized CASSCF (32,34) orbitals from our earlier work, further orbital reoptimization was not required and only Stochastic-CASCI and DMRG-CASCI calculations were performed in the localised basis.Localized orbitals have been chosen as they represent a convenient basis for analysis of electron correlation effects that involve CT excitations between the metal center and organic macrocycle.
Within the Stochastic-GASSCF method, localised active orbitals were separated into two GAS sub-spaces, namely a (18,16) subspace, consisting exclusively of (π, π * ) orbitals, and a (14,18) subspace consisting of (3d, σ N , 4s, 4p, d )orbitals.No excitations between the two active spaces were allowed (disconnected spaces), and orbitals have been relaxed under the mean-field generated by the GAS wave function [28].For the Stochastic-GASSCF procedure a convergence threshold of 1 × 10 −4 E h over the energy value was utilized.
Stochastic-GASSCF method.GAS wave functions are built by defining a number of active subspaces; within each subspace a full CI expansion is generated, while restricting the number of interspace excitations [28,42].GAS spaces are defined "disconnected " if no excitations among them are permitted, while they are defined "connected " if interspace excitations are permitted.In the same GAS wave function both connected and disconnected spaces can exist.In GASSCF, orbitals are variationally optimized under the field generated by   the GAS-CI wave function.As in the CASSCF method, all intraspace orbital rotations in GASSCF (such as GAS1-GAS1, GAS2-GAS2) are "redundant", as already described by the intra-space excitations of the GAS wave function, and can be excluded from the orbital optimization.Interspace orbital rotations (such as GAS1-GAS2), however, are only partially redundant, as inter-space excitations in the many-body expansion are absent or limited by the GAS rules.Thus, these rotations have to be considered in the orbital optimization step, in addition to inactive-active, inactive-virtual, and active-virtual rotations.Since the GAS restrictions reduce to CAS in the limit of the highest number of interspace excitations, the interspace rotations become more and more redundant for connected spaces with an increasing number of interspace excitations.This has the effect of slowing down the convergence of the GASSCF optimization.Although highly constrained, GAS wave functions with disconneced spaces are of great theoretical and practical interest.From a practical standpoint they do not suffer of the "orbital redundancies" problems discussed above.Also, one-body density matrices of disconnected GAS spaces are block diagonal and the full diagonalization of them (natural orbitals) is equivalent to block diagonalization of separate GAS subspaces ("pseudo-natural orbitals").This implies that the diagonalization of the one-body density matrix is an invariant transformation and keeps the GAS spaces separated.For connected GAS subspaces, the off-diagonal elements between orbitals belonging to different GAS spaces do not vanish and diagonalization of the one-body density matrix becomes a non-invariant rotation, because it mixes orbitals from different GAS subspaces.This implies that orbital occupation numbers are not well defined in this case.
From a chemical standpoint disconnected spaces can be utilized to selectively allow or forbid interspace excitations that can be related to specific electron cor- (2) π2 (3) π3 (4) π4 (5) π5 (6) π6  relation mechanisms, such as charge-transfer excitations.Including or removing entire classes of excitations has an impact on relative energies providing a direct and quantitative measure of the role of the class of excitations excluded.This feature will be used in this work to analyze the π backdonation mechanism in the Fe(II) porphyrin.
From an algorithmic point of view the GAS ansatz is effectively a restriction of the Hamiltonian within the corresponding complete active space, removing matrix elements that violate the GAS conditions.Thus, GAS constrains can easily be implemented within the FCIQMC framework.Only the sampling of matrix elements under GAS restriction with disconnected spaces will be described here.For a detailed description of the FCIQMC algorithm, we refer to the literature [26,27,43,44].
In FCIQMC, the sampling of the Hamiltonian reduces to randomly picking pairs of determinants |i , |j with H ij = 0 for given |i .This so-called excitation generation in practice requires a random selection of one/two electron/s and hole/s from |i for single/double excitations [26].This random sampling has to be restricted according to the GAS constraints.
For disconnected spaces where the particle number per space is constant, the excitation generator can easily and efficiently be implemented by random selection of electron(s) and subsequently of orbital(s) within the same GAS subspace.
The first step for the sampling of matrix elements is the random choice between a single or a double excitation, which is unaffected by the GAS approach.For single excitations, the GAS restrictions require the hole to come from the same GAS subspace as the electron; therefore, first an electron is chosen uniformly from all available electrons and then a hole is chosen from the same active subspace with probability, p(a|i), proportional to the matrix element H ij corresponding to the single excitation, namely where |D a i = Êai |D , represents a single excitation from the determinant, |D , and the denominator is a normalization factor with respect to all the possible single-excitations that can be generated from |D , such that all probability sum up to 1.
For double excitations, a pair of electrons is chosen uniformly, which also implies the selection of two (or one) active subspaces.From the first active subspace, a hole is then chosen uniformly.In the case of mixed-spin excitations (one α-electron and one β-electron), the spin of the hole does not have to match that of the electron in the same active space and is chosen randomly.This criterion allows for spin-flip excitations.In the Stochastic-GAS code spin-flip excitations can also be avoided forcing the hole in one GAS subspace to have the same spin projection as the electron in the same subspace.The second hole is then picked from the active space of the second electron, again with probability proportional to the corresponding double excitation matrix element The spin is chosen to ensure conservation of the total spin projection (the one defined by the reference determinant).In the case where both electrons are coming from the same GAS subspace, this scheme reduces to a heat-bath sampling [45] with the first three steps replaced by a uniform selection.FCIQMC setup.The FCIQMC setup was similar to earlier work [5] with the difference that the recently implemented adaptive shift approach [46] and in case of Stochastic CASCI, the precomputed heat bath excitation generator (PCHB), a variant of the heat bath sampling [45], were employed to speed up calculations.The adaptive shift varies the shift per non-initiator walker to cure the initiator approximation error and speeds up the convergence with respect to walker population.As the expression of the eigenvalue problem within the adaptive shift procedure is modified, also the reduced one-and two-body density matrices are to be sampled in a different way as compared to the conventional FCIQMC algorithm.For this work we used the following non-variational expression for the two-body density matrix elements where I denotes the set of all initiator determinants and Ψ total the instantaneous stochastically sampled wave function.The PCHB excitation generator samples from non-uniform probability distributions in a different manner and with different scaling characteristics than the Cauchy Schwartz excitation generator used before.[5] Unlike the adaptive shift approach, PCHB excitation generation does not introduce further approximations to the dynamics [47] and up to stochastic noise the final result for a given number of walkers is not affected, but the time per iteration is reduced.
In the CASCI (32,34) calculations the walker population was increased up to 2 × 10 9 walkers, to reach convergence, while in the Stochastic-GASSCF[ (18,16), (14,18)] calculations a much smaller population of up to 2 × 10 8 walkers was sufficient, due to the considerably smaller GASCI expansion that follows from the constrains on the CT excitations.This speedup can be estimated by comparing the different sizes of the configuration spaces.The number of determinants for a given number of electrons N , orbitals n, and spin projection S z in an active space is given by Considering that only the total spin projection is specified, and that spin projection fluctuations have been permitted within each GAS subspace, the ratio between CAS and GAS configurations is for the triplet promptly obtained as: Ω(CAS) Ω(GAS) = Ω(32, 34, 1) DMRG setup.To obtain the most accurate estimate of the exact full CI energy within the DMRG approach, a maximum bond dimension (m) value of 10000 was chosen for the natural orbital basis in D 2h point group symmetry.For the localised basis (C 1 point group symmetry) a lower m value of 4000 had to be chosen due to a larger computational cost.Each DMRG-CASCI calculation has been run with up to 20 sweeps, however the simulation was stopped after the energy difference from the previous sweep has reached 1 × 10 −6 E h .
To avoid convergence to local minima, orbital ordering according to the Fiedler vector of the mutual information matrix [48] and the perturbative correction according to Ref. 49 were employed, the latter one only for the first six sweeps in each calculation.Orbital entanglement measures, i. e. the singleorbital entropy and mutual information (see below) based on DMRG-CI calculations have been computed both for the natural and the localised orbital basis.
Orbital entanglement.Consider partitioning the system into two subsystems, |a and |b , where |a consists only of one orbital and |b of the remainder of the active space.We may obtain the one-orbital reduced density matrix (1o-RDM) by tracing out the states belonging to the environment: where c ab are the CI coefficients for configurations with a given occupation of the subsystems a and b.Since each orbital may have four possible occupation states: empty, spin-up, spin-down or doubly-occupied, ρ has a dimension of 4 × 4. We may define the single-orbital entropy s(1) i [9,10,13,14] for each orbital i as a von Neumann entropy term: where ω α,i are the eigenvalues of the 1o-RDM.The single-orbital entropy quantifies the entanglement of the orbital i with all other orbitals.Analogously, we may define for an orbital pair i, j a two-orbital RDM (2o-RDM) and the two-orbital entropy s(2) ij by placing two orbitals into the subsystem a in Eq. 6. s(2) ij describes how two orbitals are entangled with all other orbitals.From the expressions for the single-orbital entropy and two-orbital entropy, we may define the mutual information I ij as which describes the entanglement of two orbitals i and j with each other.Software.The Stochastic-CASSCF and Stochastic-GASSCF calculations have been performed within the OpenMolcas chemistry software package [50] interfaced to the NECI program [47] responsible for the FCIQMC CI optimizations.
The DMRG-CASCI calculations were performed with the QCMaquis [51,52] program, using the same molecular orbital electron integral file (FCIDUMP) produced by OpenMolcas.The entanglement diagrams were produced by the AutoCAS [53,54] program.

Results
Energies and Occupation numbers.Stochastic-CASSCF and DMRG-CASCI total energies for the 3 E g and 5 A 1g spin states, in the localised and delocalised molecular orbital basis, and Stochastic-GASSCF in the localized orbital basis are summarised in Table 1.
Table 1: Energies of the 3 E g and 5 A 1g states in Hartree and the spin gap ∆E in kcal mol −1 for the (32,34) active space in the localised and delocalised basis.The largest walker numbers for Stochastic CASSCF were 1 × 10 9 (1B) in the natural orbital basis and 2 × 10 9 (2B) in the localised one.The largest walker number for Stochastic-GASSCF[ (18,16), (14,18)] was 2 × 10 8 (200M) in the localised orbital basis.The errors in parentheses are given by the 2σ deviation of the energy calculated from the reduced density matrices.In the case of the localised basis equation 3 was used to sample the density matrices, while in the delocalised basis the contraction in equation 3 extended over all determinants.DMRG used an m value of 10 000 in the natural orbital basis and 4000 in the localised one.

State
The occupation numbers obtained from our DMRG-CASCI calculations in the delocalised basis, are in very good agreement with our Stochastic CASSCF calculation from previous work [5] and can be found in the Supporting information.
CI coefficients.The leading 10 4 Slater determinants obtained from the FCIQMC dynamic with 2 × 10 9 walkers in the localized orbital basis, were grouped into chemically relevant subsets according to the type of excitations they represent with respect to the Hartree-Fock reference determinant (see Ta-ble 2).For a set of orbitals O and P the expression (O P ) denotes the set of all Slater determinants that are generated from the reference by at least one annihilation in O and one creation in P .These subsets can be combined among each other to define new sets by forming unions, denoted with ∪, intersections, denoted with ∩.A subset relationship will be denoted with ⊂.The variable Ω denotes the set of all active orbitals in our system.It has to be emphasized, that (O P ) may contain simultaneous excitations from and into orbitals which are neither contained in O nor P .If we require that only electrons in O are annihilated and only in P created, we write pure(O P ) for the set of pure excitations.
In Table 2 the weights of these subsets have been reported as normalized values with respect to the entire wave function (p), and in an intermediatenormalization (p), with respect to the most populated reference determinant.The latter has been utilized to allow a more direct comparison of equivalent excitations in the wave functions of the two spin states considered.In fact, if two multiconfigurational wave functions with different weights on the reference determinants are to be compared, a direct comparison of the normalized weights might be misleading.Even if the weight of an excited determinant is the same in both wave functions, it represents a larger correction for that wave function with a smaller weight in the reference determinant.Intermediate normalization obviates this problem.
Orbital entanglement.Figure 4 and Figure 5 show the DMRG-CASCI orbital entanglement diagrams for the 3 E g and 5 A 1g spin states in the delocalized (natural orbitals, Figure 1 and Figure 2) and localised orbital basis (Fig. 3), respectively.
GASSCF density difference.Total electron density differences between Stochastic-CASSCF and Stochastic-GASSCF wave functions for the 3 E g and 5 A 1g spin states are displayed in Figure 6.For GAS wave functions with disconnected spaces it is possible to perform invariant orbital transformation of the GASSCF orbitals into natural orbitals (equivalent to pseudo-natural orbitals) and the electron density is then given by where the sum extends over the occupied orbitals (inactive and active), φ i refer to the natural orbitals, and ω i to their occupation numbers.For inactive orbitals ω i = 2, and for active orbitals ω i ∈ [0, 2].In the case of a GAS wave function with connected spaces, Equation 9does not hold and should be replaced with the more general form Computational Cost.Electron integral evaluation and their transformation to MO basis and Super-CI optimization of orbitals within the Stochastic-GASSCF procedure have been performed on a typical desktop workstation.The     Plot legend as in Fig. 4. The orbital numbering is the one given in Fig. 3.The two GAS subspaces used in the Stochastic-GASSCF[ (18,16), (14,18)] are highlighted with the blue and green background, respectively.

Discussion
CASCI energies are invariant under unitary orbital rotations in the active space.In the converged limit, this is also true for FCIQMC and DMRG energy estimates.Upon localisation of the Stochastic-CASSCF natural orbitals the number of determinants in the corresponding CI expansion increases as compared to the CI expansion in the natural orbital basis.This is best illustrated by considering the σ bonding and σ * antibonding molecular orbitals in the delocalised basis (Figure 1.11 and 1.12).Upon localisation into the corresponding σ 4N and 3d x 2 −y 2 orbitals (Figure 3.18 and 3.12), the atomic orbital mixing has to be described at the level of the CI expansion.Additionally, the D2h point group symmetry cannot be utilized in the localized basis; instead the C 1 point group symmetry has been employed, that leads to a further increase in the size of the CI expansion.The increased multiconfigurational character of the wave function leads to slower convergence of both FCIQMC and DMRG calculations with respect to the number of walkers and bond dimension value, respectively.
The Stochastic-CASCI calculations in the localised basis, using the adaptive shift procedure and 2 × 10 9 walkers, resulted in similar total energies and spin gaps as compared to published values in the natural orbital basis [8], where no adaptive shift technique was employed and 1 × 10 9 walkers used.A difference of 6 × 10 −4 E h between localised and delocalised basis is obtained (see Table 1), indicating the good convergence of the FCIQMC dynamics and that a higher walker population is generally necessary when covalent bondings are described by the CI expansion instead of orbital mixing.
For the DMRG calculations, the difference between energies in the localised and delocalised orbital basis is in the order of 7 × 10 −4 E h , and can be explained by considering that a smaller m value had to be employed for the localized basis, due to the increased computational requirements linked to the C 1 point group symmetry utilized.
The loss in accuracy under orbital localisation is similar for the FCIQMC and DMRG procedures, and they yield the same spin gaps in the same orbital basis (see Table 1).Not only the energy, but other properties of the wave function, such as the natural orbital occupation numbers, are very similar for the two methods.
Although FCIQMC a DMRG algorithms are conceptually different, they are qualitatively and quantitatively in agreement in the description of the Fe(II) porphyrin model system within the (32,34) active space.This agreement allows to compare FCIQMC and DMRG properties directly, and analyze orbital entanglement and CI coefficients of various classes of excitations on equal footing.
Natural Orbital Basis.The most prominent feature in both 3 E g and 5 A 1g states are the Gouterman's orbitals 16, 24, 29, and 33 (see Fig. 1), which show the strongest entanglement for both spin states (Figure 4).As already anticipated by Gouterman [55], these orbitals play an important role within the π system, and as previously discussed by us [5], they are crucial in differentially stabilizing the triplet over the quintet spin states.The entanglement diagrams effectively confirm the crucial role of the Gouterman orbitals.However, the entanglement diagrams also show that other π orbitals are strongly correlated (e. g. orbitals 21, 27 and 34) as they are correlated to one or more Gouterman's orbitals and show large values of single-orbital entropy.Generally speaking the entire system of π orbitals is important to capture the important correlation effects related to the aromatic macrocycle.
From the entanglement diagrams of Figure 4 another aspect of the triplet spin state stabilization, already discussed in Ref. 5, is highlighted, namely the breathing effect of the 3d electrons.This is exemplified by the entanglement between orbitals 22, 23, 24 and 26, 29, 30, respectively.At the first glance, the entanglement within these orbital triples is similar in both states and not very strong.However, orbital 26, which is relevant for the orbital breathing phenomenon, along with orbitals 27, 28, and 31 in the natural orbital basis, are localised in the 5 A 1g state but delocalised in the 3 E g spin state (see Figure 2).Thus, it is not possible to directly compare the entanglement diagrams of the 3 E g and 5 A 1g spin states and draw conclusions on the role of orbital breathing and charge transfer processes.The role of orbital breathing and CT mechanisms are better illustrated by the entanglement diagrams in the localised orbital basis (Figure 5), discussed in the next section.
Localised Orbital Basis.As shown in Table 2, the normalized weight (p) of the reference determinant of the triplet is substantially smaller than the normalized weight of the reference determinant of the quintet spin state, indicating that correlation effects are in general more pronounced for the triplet.We will discuss metal centered correlation effects first and CT effects between metal and π system later.
The σ 4N and 3d x 2 −y 2 orbitals (12 and 18 in Figure 3).have the largest entanglement (see Figure 5), and describe the σ bonding interaction between the metal center and the macrocycle.In the quintet state the σ-bond entanglement is less pronounced, providing a simple route for an explanation of the geometrical differences between triplet and quintet spin states, with the quintet states in general characterized by larger Fe-N bond distances.This result identically emerges from the inspection of the CI coefficients listed in Table 2.In fact, considering the normalized weights (p), the set of excitations from σ N to 3d orbitals makes up 33 % of the total wave function in the triplet, and 12 % in the quintet state and are dominated by the (σ 4N 3d x 2 −y 2 ) excitations.As expected, the correlation within the metal 3d orbitals (orbitals 14-18 in Figure 3) is stronger in the triplet spin state, and it is clearly shown both by the orbital entanglement diagrams and the CI weight of (3d 3d) excitations.Both the 3d correlation and σ donation are already contained in the (8, 6) minimal active space.The double shell correlation of 3d orbitals with the correlating d functions (19 to 23 in Figure 3) accounts for the 2.72 % and 1.40 % weight for the triplet and quintet state, respectively, and represents an important factor for a balanced description of the wave functions of the Fe(II) porphyrin, that is captured in the (8,11) active space utilized by Pierloot and Co-workers.[56] Charge-transfer excitations between metal orbitals and π orbitals of the macrocycle are only contained in the large (32,34) active space.These excitations differentially favour the triplet over the quintet.The entanglement between the π 9 and 4p z orbitals (9 and 27 in Figure 5) is the second largest entanglement value with a corresponding high CI weight for (π 4p z ) in both spin states.Although the normalized weight is slightly higher for the quintet spin state, with 19.40 % and 22.87 %, for the weight of the triplet and quintet, respectively, in the intermediate normalization this excitation has a larger im-pact on the wave function of the triplet spin state.The 4p z orbital is empty in the reference determinant and has the correct symmetry to be used as correlating orbital for the electrons of the π system.It is to be expected that this excitation will reduce in weight if correlating orbitals were included into the active space.This aspect is currently under investigation.
As the entanglement diagrams (Figure 5) and the CI weight of (3d π * ) show (Table 2), there is significant π backdonation from metal 3d orbitals into the π system.This correlation is much more pronounced in the triplet, hence having a major impact on the spin gap.The π backdonation mainly involves the Gouterman's orbitals which is illustrated by the CI weight of the restricted subset (3d π * G ) ⊂ (3d π * ).This special role of the Gouterman's orbitals is also manifested in their single orbital entropy, which is larger than of the rest of the π orbitals.On the other hand it has to be emphasized, that the CI weight of the pure excitations pure(3d π * G ) is significantly lower than the case where they are coupled to other excitations, which illustrates that the π backdonation is not a pure excitation from 3d into Gouterman's orbitals, but a concerted mechanism.In our previous work [5] there was already evidence for this observation because natural orbital occupation numbers of Gouterman's orbitals were close to 0.0 and 2.0 when they were added to an active space consisting of solely metal orbitals.Only when the full π system was correlated, the occupation numbers of the Gouterman orbitals changed to approximatively 1.7 and 0.3.[5] The "breathing" of the correlated 3d electrons enables an other complex correlation mechanism that involves metal 3d-orbitals, correlating d orbitals and π * orbitals.This mechanism is also captured by the entanglement diagrams (Figure 5) and corresponds to the entaglement of the 3d  5) shows that the orbitals are more entangled in the triplet state.Also in the CI weights (Table 2) the "breathing" manifests itself in the weight of the excitations in (3d While the weight of this subset is similar for both spin states, the intermediatenormalized weight of 7 % shows, that the "breathing" is only for the triplet a substantial correction to the reference wave function. The π → d ligand-to-metal charge-transfer excitations are important contributions to the overall wave function as illustrated by the CI weight of (π d ) (Table 2).Figure 5 shows entangled orbital pairs π 6 − d xz (6 − 20) and π 7 − d yz (7 − 21) as most relevant for these excitations, which provide a mechanism for the charge delocalisation similar to orbital breathing but working in the opposite direction, a mechanism that reduces correlation at the level of the macrocycle via delocalization of charge.
Stochastic-GASSCF In the previous paragraph we have discussed the contribution of correlation effects that involve metal and π orbitals by inspection of orbital entanglement and CI weights of various classes of excitations.In this sec-tion the Stochastic-GASSCF approach will be utilized to quantitively determine the impact of such mechanisms on spin-state energetics.
As it is expected from the Hylleraas-Undheim-MacDonald theorem, the variational energy is increased upon reduction of the configuration space that follows the GAS restrictions (see Table 1).[57,58] More interesting is the fact, that the spin gap is reduced by 2.5 kcal mol −1 showing that the triplet is stabilized by correlation effects among electrons of the metal center and of the π system.These effects are similar in magnitude to the semicore correlation effects reported by Phung et al. [7] and Li Manni et al. [8].
However, the fact that the Stochastic-GASSCF[ (18,16), (14,18)] still favors the triplet spin state, suggests that the stabilization of the triplet is not exclusively a consequence of many-body correlation effects; it also relates to meanfield effects of the correlated macrocycle.In fact, in the GASSCF procedure, although interspace excitations have been forbidden, the intra-space excitation within the π-system are still considered, and they are capable to change the mean-(ligand-)field effect onto the metal center.This is to be compared to the case where the π system is left outside the active space and treated at the single configurational level of theory.This justifies why, despite CT configurations are removed, the triplet is still lower in energy.
The impact of GASSCF constraints can also be observed in the electron density difference between the CASSCF and GASSCF wave functions ∆ρ 1 := ρ(CASSCF) − ρ(GASSCF) (see Figure 6).In going from GASSCF, where CT configurations are excluded from the wave function, to CASSCF, where all possible configurations are allowed, the density increases in the π system, and decreases in the metal.This in turn polarizes the nitrogen and enhances sigma donation.These effects are present in both spin states, but more pronounced in the triplet.A similar analysis was performed in our previous work [8] where the density difference between the Stochastic-CASSCF (32,34) and the restricted open-shell Hartree-Fock (ROHF) wave function ∆ρ 2 := ρ(CASSCF)−ρ(ROHF) was considered.As qualitatively expected, the overall difference between the uncorrelated ROHF and CASSCF density is higher than between GASSCF and CASSCF ( |∆ρ 2 (r)| dr > |∆ρ 1 (r)| dr).The polarization of the nitrogen relative to the density change in the π sytem is larger in the case of ∆ρ 1 , while the nitrogen polarization and density changes in the π system are of similar magnitude in the case of ∆ρ 2 .This can be explained by the fact that (π π * ) excitations are already contained in the GASSCF, but not in the uncorrelated ROHF wave function.This intraspace correlation affects the mean field of the macrocycle onto the metal center and differentially stabilizes the triplet as discussed above.

Conclusion
In this work, Stochastic-CASSCF, Stochastic-GASSCF and DMRG-CASCI procedures have been applied to an Fe(II) porphyrin model complex in a localized orbital basis, to understand the correlation effects that govern the spin gap between the 5 A 1g and 3 E g electronic states.The choice of the localized orbitals has been motivated by the possibility to conveniently separate local excitations from CT excitations between metal center and macrocycle.
The FCIQMC and DMRG approaches have provided insights that are in qualitatively and quantitatively good agreement, with a difference of 5 × 10 −4 E h in total energies, a negligible difference of < 1 × 10 −1 kcal mol −1 in the spin gap estimates, and very similar natural orbital occupation numbers.
From the two methods, two complementary pieces of information on the role of electron correlation effects in metal-porphyrins were obtained.With FCIQMC it has been possible to directly analyse the leading electronic configurations that characterize the multi-configurational wave function of the 3 E g and 5 A 1g spin-states, while the DMRG procedure, in the same orbital basis, has produced entanglement diagrams that facilitate the understanding of complex electron correlation effects in metal-porphyrins.
It has been shown that multi-orbital correlation mechanisms between the metal and π system, such as the 3d d π * "breathing" processes are required for a qualitatively correct description of the wave function, and they differentially stabilize the intermediate spin state over the high spin states.
A remarkable finding is about the role of the Gouterman's orbitals.It has been demonstrated that the π backdonation 3d π * happens mainly into the Gouterman's orbitals (large values of orbital entanglement), but requires simultaneous excitations in the whole π-system.This feature is at the core of the necessity to include the whole π system of the macrocycle into the active space, and at the core of the failure of second order perturbation procedures described in earlier work.
The Stochastic-GASSCF method with disconnected spaces has been developed within the FCIQMC framework, to quantify the impact of charge-transfer excitations.Within this method, localized metal centered orbitals and π orbitals of the macrocycle were separated into two active subspaces.With no interspace excitations, this setup effectively removed the stabilizing CT configurations from the many-body expansion of the two spin states considered, responsible for the π backdonation mechanism, and reduced the spin gap from −3.5 kcal mol −1 to −1.0 kcal mol −1 .
The density difference between the Stochastic-CASSCF and Stochastic-GASSCF wave functions shows that CT correlation effects between the metal and macrocycle amplify 3d π * backdonation and σ donation of nitrogen ligands, and in doing so, confirming the suggested correlation enhanced σ donation / π backdonation mechanism.
One important question is whether the results obtained in this work can be transferred to the real Fe(II) porphyrin system, with its pyrrole β-carbon atoms, and spin gaps predictions in adiabatic conditions.The model system without β-carbons utilized in this work contains a planar and fully conjugated macrocycle, that mimics the chemical environment around the metal center of real metalporphyrins.The β-carbons in Fe(II) porphyrin, increase ring-correlation at the level of the π-system of the macrocycle.Thus, further stabilization of the triplet over the quintet is to be expected when considering systems that contain the pyrrole β-carbons.
Adiabatic conditions will lower the energy of the quintet state with respect to the same state at the triplet geometry.The correlation mechanisms described here will still play an important role in differentially stabilizing the intermediate spin state.Nevertheless it might be necessary to include even more correlation effects for a correct description in adiabatic conditions, considering that convergence of the spin-gap with respect to the active space size has not been reached within our calculations, as shown in Figure 2 of Reference [5].The identification and quantification of further correlation effects are currently being investigated.
Future work on Stochastic-GASSCF will develop in two directions.One goal is to enable interspace excitations which can be used in a similar but more flexible way as the disconnected spaces in this work.Another goal is to apply the approximate, precomputed heat-bath excitation generation and other more efficient sampling strategies in the GAS framework for enhanced performance [45,59].Considering the reduced computational costs associated with Stochastic-GASSCF as compared to the corresponding Stochastic-CASSCF wave function optimization, larger active spaces are within reach, providing practical routes for investigating correlation mechanisms.

Table of Contents text
The role of charge transfer electron correlation mechanisms that are responsible for the differential stabilization of the 3 E g over the 5 A 1g state in a Fe porphyrin model system, have been studied by Stochastic-CASSCF, DMRG and a preliminary version of the Stochastic-GASSCF method.Orbital entanglement diagrams and CI coefficients of the wave function in a localised orbital basis allow for an effective interpretation of charge transfer configurations.Within Stochastic-GASSCF, excitations between metal and ligand orbitals have been selectively removed from the complete CI expansion.These excitations are key to the differential stabilization of the triplet and once removed the spin-gap noticeably reduced.

Figure 2 :
Figure 2: Natural active orbitals for the 3 E g state of the Fe(II) porphyrin model system.Only the orbitals that differ from the 5 A 1g state orbitals from Figure 1 are depicted.

Figure 4 :
Figure 4: DMRG-CASCI entanglement diagrams for the (a) 3 E g and (b) 5 A 1g spin states, in the natural orbital basis of the Fe(II) porphyrin model.The area of the grey circles is proportional to the single-orbital entropy s(1) i , the lines connecting the orbitals show their mutual information value I ij .Solid lines denote I ij values of at least 0.1, dashed lines I ij values of at least 0.01.The orbital numbering is the one used in Fig. 1 and 2.

Figure 5 :
Figure 5: DMRG-CASCI orbital entanglement diagrams for the (a) 3 E g and (b) 5 A 1g spin states of the Fe(II) porphyrin model in the localised orbital basis.Plot legend as in Fig.4.The orbital numbering is the one given in Fig.3.The two GAS subspaces used in the Stochastic-GASSCF[(18,16),(14,18)] are highlighted with the blue and green background, respectively.
Figure 1: Natural active orbitals for the 5 A 1g state of the Fe(II) porphyrin model system sorted by irreducible representation and occupation number.The four Gouterman π-frontier orbitals are labelled with the additional subscript G.

Table 2 :
[5]cent weight, p, of different subsets of the 10 4 most populated Slater determinants in the CI expansion for the localised basis.The intermediately normalized weight p = p/p ref is normalized onto the weight of the reference determinant.The first four sets of excitations do not involve π-orbitals of the macrocycle, and are already contained in the smaller active spaces discussed in Reference[5].The subsequent sets involve π-orbitals of the macrocycle.The variable π G indicates the four Gouterman π orbitals, π 8 , π 9 , π * 1 , π * 2 (Figure3.8,3.9, 3.28, and 3.29).If a set of excitations is a subset of another set in the table, it is indented.