Spin-pure Stochastic-CASSCF via GUGA-FCIQMC applied to Iron Sulfur Clusters

In this work we demonstrate how to compute the one- and two-body reduced density matrices within the spin-adapted full configuration interaction quantum Monte Carlo (FCIQMC) method, which is based on the graphical unitary group approach (GUGA). This allows us to use GUGA-FCIQMC as a spin-pure configuration interaction (CI) eigensolver within the complete active space self-consistent field (CASSCF) procedure, and hence to stochastically treat active spaces far larger than conventional CI solvers whilst variationally relaxing orbitals for specific spin-pure states. We apply the method to investigate the spin-ladder in iron-sulfur dimer and tetramer model systems. We demonstrate the importance of the orbital relaxation by comparing the Heisenberg model magnetic coupling parameters from the CASSCF procedure to those from a CI-only procedure based on restricted open-shell Hartree-Fock orbitals. We show that orbital relaxation differentially stabilizes the lower spin states, thus enlarging the coupling parameters with respect to the values predicted by ignoring orbital relaxation effects. Moreover, we find that while CI eigenvalues are well fit by a simple bilinear Heisenberg Hamiltonian, the CASSCF eigenvalues exhibits deviations that necessitate the inclusion of biquadratic terms in the model Hamiltonian.


I. INTRODUCTION
The complete active space self-consistent field (CASSCF) method is a well-established approach in quantum chemistry for the treatment of strongly correlated electron systems with substantial multi-reference character. [1][2][3][4][5][6][7][8] Important static correlation effects are rigorously described within the active space, consisting of the most important orbitals and electrons, while the effect of the environment (electrons not included in the active space) is accounted for at the mean-field level via a variational orbital optimization (the SCF procedure). One-and two-body reduced density matrices (1-and 2-RDMs) within the active space are necessary to perform orbital rotations between the active orbitals and the environment, whether a secondorder Newton-Raphson formulation, 2,7,[9][10][11] or the simplified Super-CI technique with an average Fock operator is utilized. 4 If applicable, exact diagonalization techniques [12][13][14] are utilized to obtain eigenvalues, eigenvectors and the RDMs associated to the CAS configuration interaction (CASCI) Hamiltonian. However, due to the exponential scaling of CASCI with respect to the size of the active space, exact diagonalization techniques are restricted to at most about 18 electrons in 18 orbitals, CAS(18e,18o), on serial architecture. 15,16 More recent massively parallel implementations allow sizes up to CAS(24e,24o). 17 Another strategy is to use methods that approximate the full-CI wave function in the active space, like the density matrix renormalization group approach (DMRG), [18][19][20][21][22][23][24][25][26][27][28][29][30] full configuration interaction quantum Monte Carlo (FCIQMC), 6,31-35 selected configuration interaction (Selected-CI) approaches, [36][37][38][39][40][41][42][43][44][45] recently implemented in a spin-adapted form [46][47][48] , as a) Electronic mail: w.dobrautz@fkf.mpg.de well as the generalized active space approach, [49][50][51] as CI-eigensolvers within the CASSCF framework. These approaches allow the study of much larger active spaces. 6,39,[51][52][53][54][55][56][57][58][59] The use of FCIQMC as the CASSCF CI-eigensolver within the Super-CI framework, termed Stochastic-CASSCF, 6 has been developed in our group, and used to study a number of strongly correlated systems, such as model systems of Fe(II)-porphyrins and the correlation mechanisms that differentially stabilize the intermediate spin states over the high spin states, 51,52,54 and model systems of corner-sharing cuprates. 53 The original Stochastic-CASSCF implementation was formulated using Slater-Determinants (SDs) as many-body basis for the CASCI wave function expansion. As SDs are not necessarily eigenfunctions of the total spin operator, its applicability is bound to the intrinsic spin structure of the system studied. If the lowspin states are energetically more stable and well separated from higher spin states, it is possible to obtain essentially spin-pure wave functions when using an SD basis. However, when high-spin states are more stable than low-spin states, and/or a number of spin states are nearly degenerate, it is very difficult to obtain spin-pure solutions, or target states other than the ground state, with a SD basis.
In this paper we present an algorithm for the calculation of 1-and 2-RDMs within the spin-adapted implementation of FCIQMC via the graphical unitary group approach (GUGA-FCIQMC). 60 GUGA-FCIQMC has been implemented within the NECI code, 31,61 and provides accurate spin-adapted wave functions and RDMs for active space sizes out-of-reach for conventional exact CI-eigensolvers. 57,58 As already done for the original Stochastic-CASSCF 6 , the sampled 1-and 2-RDMs are then utilized within the Super-CI procedure as implemented in the OpenMolcas chemistry software package, 16 to perform the orbital relaxation step. Thus, via the interface of the NECI code and OpenMolcas 16 , it is possible to perform spin-adapted state-specific (or state-average, if RDMs of different states are weightedaveraged prior the Super-CI step) Stochastic-CASSCF optimizations, targeting any desired spin state. The spin-pure Stochastic-CASSCF allows us to obtain variationally optimized molecular orbitals, which in turn enable the calculation of spin gaps, unbiased from the choice of the starting orbitals.
The applicability and the importance of the method is shown through the investigation of the spin ladder of two iron-sulfur (FeS) clusters. Poly-nuclear transitionmetal (PNTM) clusters are of major importance in organometallic chemistry and as cofactors in biology, and are involved in a multitude of processes, including photosynthesis, respiration and nitrogen fixation [62][63][64] , being responsible for redox reactions [65][66][67] and electron transfer, [68][69][70][71][72][73][74][75][76] act as catalytic agents and even provide a redox sensory function. 77 A theoretical understanding of the intricate interplay of the energetically low-lying spin states of these systems, guided by accurate numerical results, could provide insights towards the synthetic realization of these processes. Especially, because direct experimental measurements targeting the electronic structures of these systems are often hindered by the large number of overlapping electronic states, and corresponding vibrational modes at finite temperatures. [78][79][80] In addition, some energetically low-lying excited states are inaccessible by accurate optical absorption experiments due to being electric-dipole forbidden transitions. 81 Spin-pure stochastic RDM sampling allows us to formulate a spin-adapted Stochastic-CASSCF and gives us access to properties encoded in the 1-and 2-RDMs, such as spin-spin correlation functions. Using CASSCF wave functions of various active space size and composition we will study and discuss how spin gaps are affected by orbital relaxation effects. Additionally, the ab initio energies will be mapped to the (biquadratic) Heisenberg spin model [82][83][84][85][86][87][88][89][90][91][92] to show the effect of active space size and orbital relaxation on the extracted magnetic coupling parameters, which are in turn compared to the available experimental data 93,94 and other computational studies. 95,96 The remainder of this paper is organized as follows: In Sec. II we summarize the spin-adapted GUGA-FCIQMC method and in Sec. III we describe the sampling algorithm of spin-free RDMs. In Sec. IV we discuss ab initio CASSCF spin gaps and spin-spin correlation functions for an iron-sulfur dimer, Fe 2 S 2 , 57 for different active space sizes and starting orbitals, and for an [Fe 4 S 4 ] tetramer model system. We also map our ab initio results to a (biquadratic) Heisenberg model Hamiltonian, discuss the role of the CASSCF procedure when extracting the exchange parameter(s), and compare the magnetic coupling constants extracted from our computations to experimental and theoretical references. Finally, in Section V we summarize our findings and offer a general discussion on the presented topic.
An appendix is available, where we derive necessary formulas for local spin measurements (Appendix A), and spin-spin correlation functions from RDMs (App. C). We additionally supply coordinate and orbital files, compu-tational details, and comparisons with available exact results for small active spaces, a table with the data used in Fig. 10, a study on improved convergence due to stochastic noise, the protocol how we compared the orbitals in Fig. 11, details on interface and the RDM storage convention in OpenMolcas and a quick access literature overview of computation results for the Fe 2 S 2 system in the supporting information (SI). 97

II. GUGA-FCIQMC
In this section we briefly summarize the main details of the GUGA-FCIQMC implementation. More theoretical and technical aspects of the algorithm are available in the literature. 60,98 The spin-adapted implementation of the FCIQMC algorithm relies on the unitary group approach (UGA), 99,100 pioneered by Paldus, and its graphical extension (GUGA), introduced by Shavitt. 101,102 GUGA provides an efficient-to-use spin-adapted many-body basis, based on the spin-free formulation of quantum chemistry. 103 The spin-free form of the electronic Hamiltonian is given byĤ with the spin-free excitation operators,Ê ij = σâ † iσâ jσ andê ij,kl =Ê ijÊkl − δ jkÊil , defined in terms of the creation and annihilation operatorsâ † iσ ,â jσ with spatial orbitals i, j and spin σ. t ij and V ijkl represent the oneand two-electron integrals in a molecular orbital basis, and n indicates the total number of spatial orbitals.
The name unitary group approach comes from the fact that the operatorsÊ ij fulfill the same commutation relations as the generators of the unitary group of order n, U (n). 99 Paldus 104,105 identified a very efficient construction of a spin-adapted basis tailored for the electronic structure problem, based on the Gel'fand-Tsetlin basis, 106-108 a general basis for any unitary group U (n). The graphical extension, GUGA, provides an efficient way to calculate Hamiltonian matrix elements, ν|Ĥ|µ between different CSFs, |µ and |ν , within a chosen spin-symmetry sector. The combination of an efficient protocol for computing Hamiltonian matrix elements, and storage of the CI coefficients, enables an effective spin-adapted formulation of exact CI eigensolvers, such as CAS 1 and GAS 49 , perturbation theory methodologies, such as CASPT2 109 , GASPT2 110 and SplitGAS 111 as well as the FCIQMC approach within the GUGA framework. 60 The FCIQMC algorithm 33,34 is based on the imaginary-time (τ = it) Schrödinger equation, (2) which, after formal integration and a first-order Taylor expansion yields an iterable expression for the eigen-state, |Ψ(τ ) FCIQMC stochastically samples the FCI wave function, |Ψ(τ ) , of a system by a set of so-called walkers and yields estimates for the ground-and excitedstate 112 energies and properties 113 via the one-and twobody RDMs. 114 Theoretical and algorithmic details on FCIQMC can be found in the literature, 33,34 especially in the recently published review article [31]. At the heart of the FCIQMC algorithm is the so called spawning step, which stochastically samples the off-diagonal contribution to the imaginary-time evolution of the targeted state, with c ν (τ ) being the coefficient of basis state function |ν at the imaginary-time τ , of the FCI expansion |Ψ(τ ) = ν c ν (τ ) |ν and p gen (µ|ν) is the so-called generation probability of choosing configuration |µ given |ν . During a FCIQMC simulation only coefficients that are at least occupied by a chosen minimum number of walkers (usually set to be the real number 1.) are kept in memory. The off-diagonal contribution in Eq. (4) is then approximated by allowing each walker on each occupied configuration |ν to spawn new walkers on configuration |µ with a non-zero Hamiltonian matrix element µ|Ĥ|ν . The process of suggesting a new configuration |µ given |ν , called the excitation generation step, is of utmost importance.
The maximal usable time-step of the simulation is limited by the relation to ensure stable dynamics. Hence, for large |H µν |/p gen (µ|ν) ratios, the time-step of the calculation, ∆τ , has to be lowered to ensure a stable simulation and is the motivation for optimizing the excitation generation step. Several schemes to obtain a close-to-optimal balance of computational effort and matrix element relation have been developed, see References [31,51,[115][116][117]. The spawning step is schematically shown in Fig. 1.
For reasons of interpretability, control, and improved convergence properties, a spin-adapted implementation of FCIQMC was long-sought after. 118 GUGA allows an efficient spin-adapted FCIQMC implementation, by constructing spin-symmetry allowed excitations as stochastic walks on the graphical representation of CSFs, the so-called Shavitt graph, as depicted in Fig. 2, and explained in more detail in References [60,98].
The GUGA allows both an efficient on-the-fly matrix element calculation and a way to select excitations from CSF |µ → |ν , and ensures the approximate relation p gen (ν|µ) ∝ |H µν |, via a so-called branching tree approach. The stochastic GUGA excitation process for a  The loop contributing to the coupling coefficient, µ|Êij|µ , is indicated by the gray area. Following Shavitt's convention the CSFs are drawn from bottom to top. (b) Exchange excitation example, for the same CSF |µ , which shows that different index combinations for exchange excitations,êij,ji andê i j ,j i , can lead to the same transition |µ → |µ = |u, 0, d, u, u, 0 , with a non-zero coupling coefficient. single excitation,Ê ij , is schematically depicted in Fig. 3, with the CSFs drawn from top to bottom. For a given CSF, |µ , and two spatial orbitals, i and j, which are chosen with a probability weighted according to the magnitude of their integral contributions, at each open-shell orbital k within the range i → j, an allowed path is chosen randomly. This process is weighted with the socalled probabilistic weight, of the remaining decision tree below the current orbital k, which ensures the desired relation p gen (ν|µ) ∝ |H µν |. Interested readers are referred to Refs. [60,98] for more details.
This stochastic process, additionally circumvents the bottleneck given by the exponentially growing connectivity between CSFs with respect to the number of openshell orbitals. Hence, our GUGA-FCIQMC method is Schematic representation of the branching-tree approach to allow efficient on-the-fly excitation generation and matrix element calculation entirely in the space of CSFs without any reference to SDs. Example given for a single excitationÊij from a CSF |µ = |u, 0, u, d, 0, u, 0 to one other |ν = |u, u, d, u, 0, 0, 0 . In the shown example, an electron is excited from the singly occupied orbital 6 in |µ to the empty orbital 2 (indicated by the arrow on the left). Spin-allowed excitation pathways are indicated by solid lines. During the random excitation process in GUGA-FCIQMC, a spin-symmetry allowed path is chosen at random, weighted according to the resulting coupling coefficient, µ|Eij|ν (indicated by the red pathway). In general, empty starting orbitals and singly occupied orbitals in the excitation range allow for two possible spin couplings (u/d). Spin-symmetryforbidden paths are indicated by the crossed-out nodes.
able to treat systems with more than 30 open-shell orbitals and due to the spin-pure formulation, it allows to specifically target any spin-symmetry sector, removes any spin-contamination, reduces the Hilbert space size, and speeds up convergence in systems with neardegenerate spin states. 60,98 However, compared to the SD-based FCIQMC, the calculation of (spin free) RDMs in GUGA-FCIQMC is considerably more challenging, and hence such RDMs have not been available until now. This has prevented to access properties, and using it as a spin-pure CI-solver within the Stochastic-CASSCF method. 6

III. GUGA-RDMS
In this section the theoretical and algorithmic details of the stochastic sampling of RDMs within the GUGA-FCIQMC method are discussed.

A. Theoretical considerations
Unbiased RDM sampling within the FCIQMC algorithm, whether in SD or CSF basis, is made possible by the replica method, where two independent dynamics are simultaneously carried to remove a strictly positive bias due to stochastic fluctuations for the diagonal RDM contributions. 114 In a SD based implementation the 1-particle RDM entries are derived from the stochastic coefficients c A I (τ ) and c B J (τ ) of two statistically independent calculations, A and B. The two-body RDMs are obtained in a similar way. For SDs, the terms I|a † iσ a jσ |J are promptly given by the well-known Slater-Condon rules. We make use of the fact that |I → |J transitions are already performed in FCIQMC during the stochastic spawning step. Hence, we reuse the information, already required for a normal simulation, to additionally sample the 1and 2-RDM elements.
In the original SD based implementation this is done by additionally storing information of the parent SD, |J , along with the spawned new SD, |I , including the parent SD encoded in a bit representation, its coefficient, in what run (A or B) this spawn happened, and other implementation specific flags.
In a parallel high-performance computing (HPC) environment, the occupied determinants are distributed among the different processors. Hence, the newly spawned walkers are kept in an array, which has to be communicated to the corresponding processor, where the newly spawned state is stored, to update the corresponding coefficients.
The spin-free one-and two-body RDMs in terms of unitary group generators 119 are defined as (following the convention of Helgaker, Jørgensen and Olsen 3 ) withÊ † ij =Ê ji , and with i, j, k, l denoting spatial orbitals, |µ and |ν being configuration state functions (CSFs) and c µ and c ν their coefficients in the ground state wave function expansion, |Ψ .
The diagonal terms of the RDMs are accumulated explicitly, and the diagonal 1-RDM terms reduce to where n i is the occupation of the spatial orbital i, which can assume the values 0, 1, or 2. The diagonal 2-RDM elements are defined as which for i = j yields and for i = j Eqs. (11) and (12) are simply products of orbital occupation numbers and the coefficients c A/B µ from two statistically independent simulations, due to the above mentioned positive bias in diagonal RDM entries. Exchange type elements of the 2-RDM, Γ ij,ji , also have diagonal contributions from the wave function which are also sampled explicitly. The detailed form of the coupling coefficients can be found in the literature. 60,102 These exchange-like terms do, however, also have off-diagonal contributions, µ|ê ij,ji |ν , which are explained in general in the following section.

B. Off-diagonal RDM entries -computational implementation and cost
Similar to the SD based RDM sampling, for each sampled RDM element we store the parent state, |µ , its coefficient, c µ , and the replica index (A or B). However, there are some important differences in the GUGA based RDM sampling compared to a SD-based implementation: (a) The one-electron coupling coefficients µ |Ê ij |µ , and the corresponding two-body terms, µ|ê ij,ji |ν , do not follow the Slater-Condon rules as for SDs. Shavitt and Paldus derived an efficient product form of these coupling coefficients, exemplified by a single excitation as, where W is a function of the step-values, d k = {0, u, d, 2}, of spatial orbital k of the step-vector representation of the two CSFs, |µ and |µ , and the intermediate value of the total spin, S k , in the cumulative sense. The step-values, d k , encode if a spatial orbital is empty, d k = 0, positively spin-coupled ∆S k = +1/2, d k = u, negatively spin-coupled, d k = d or doubly occupied, d k = 2. CSFs can be represented graphically (see Fig. 2a) where different step-values are indicated by a different tilt of the segments, and Shavitt showed that the value of the coupling coefficients only depends on the loop shape enclosed by the two coupled CSFs. Their explicit calculation scales with the number of spatial orbital indices between i and j. However, we calculate this quantity on-the-fly, during the excitation generation step, and thus, we can reuse it with no additional computational cost in the stochastic RDM sampling.
(b) Identifying the type of excitation and the involved spatial orbitals (i, j, k, l), when coupling CSFs, is a more complex operation than for SDs. CSFs can also differ in the open-shell spin-coupling, and not only in the specific spatial orbitals (i, j, k, l), yet still have a nonzero coupling coefficient. For example, in Figure 2a we show Shavitt's graphical representation of the CSF |µ = |u, 0, u, d, u, 0 (already used in Fig. 3) as the orange solid line and an excited CSF |µ = |u, 0, 2, u, 0, , 0 as the blue dashed line. Following Shavitt's convention the CSFs are drawn from bottom to top. Only the gray loop area enclosed by both CSFs contributes to the coupling coefficient, µ|E ij |µ . The two CSFs are connected by an excitation of an electron from orbital j = 5 to i = 3, indicated by the arrow. However, as one can see in Fig. 2a, the two CSFs |µ and |µ do also differ in the spin coupling of orbital k = 4 with d k = d, while d k = u (in the step-value notation). Hence, it is not as simple as performing bit-wise logical operations on alphaand beta-strings as it is possible for SDs, 120 to identify the involved spatial indices and type of excitation. We do have optimized routines to perform this excitation identification for arbitrary CSFs in our GUGA-FCIQMC code NECI, 31,61 and similar to the above-mentioned coupling coefficients, we already have the necessary information in the excitation process, within the spawning step.
(c) Certain excitation types, such as the exchangelike excitations,ê ij,ji andê ij,jk , can have multiple nonunique spatial orbital combinations leading to the same type of excitation |µ → |µ . This stems from the fact that certain contributions to the two-body coupling coefficients, µ|ê ij,kl |µ , are non-zero for alike open-shell step-values, d o = d o , above and below the loop spawned between |µ and |µ , see Figure 2b and Shavitt. 102 For example, for a pure exchange-type excitation, e ij,ji , as depicted in Figure 2b, only the spin-coupling of the open-shell orbitals differs, but there is no change in the orbital occupation. To calculate the Hamiltonian matrix element, µ|Ĥ|µ = i =j V ijji µ|ê ij,ji |µ , one needs to consider all non-zero contribution to the coupling coefficient, from orbital i below and j above the loop. Additionally, as the specific spatial orbitals, i, j (k, l) are chosen first in the excitation generation in FCIQMC, it is necessary to also take into account the possibilities that the other contributing orbitals, p(i , j ), would have been picked (as their choice could have led to the same excitations), to assign a unique total generation probability, p gen (µ |µ). However, for a correct RDM sampling we have to retain the original probability p(µ → µ |i, j, k, l) to sample a specific Γ ij,kl entry to avoid a possible double counting. Conveniently, similar to the cases (a) and (b) mentioned above, we already have access to this specific quantity, obtained during the excitation generation process and do not need to explicitly recalculate it for the stochastic RDM sampling.
The three additional necessary quantities discussed above, namely the coupling coefficient, µ |Ê ij |µ or µ |ê ij,kl |µ , the excitation type, and the probability p(µ → µ |i, j, (k, l)) are already computed in the random excitation process. Consequently, the main change to enable spin-free RDM sampling within GUGA-FCIQMC is to communicate these three additional quantities, along with the already communicated information of the parent state, |µ , its coefficient, c µ , and the replica index, A/B.
An important algorithmic advancement and routinely used feature of FCIQMC was the semi-stochastic method, 121,122 where some chosen part of the Hilbert space -usually the N D most occupied states -is treated explicitly. This is achieved by constructing the full Hamiltonian matrix H µν , ∀µ, ν ∈ {N D } and performing the imaginary-time evolution exactly. This necessitates also a change to the RDM sampling, since the RDMs contributions from states within the semistochastic space are not covered in the random excitation process anymore. These RDM contributions are treated exactly, greatly increasing their accuracy on the one hand, but -especially in the spin-free case -also increasing the computational effort. In this case it is not possible to avoid the explicit excitation identification and coupling coefficient and original generation probability calculation in GUGA-FCIQMC. However, there is only marginal computational overhead of around 10-20% associated with the spin-free RDM sampling compared to a standard two-replica FCIQMC calculation (see the SI for details).
The spin-adapted Stochastic-CASSCF method has been made available in the OpenMolcas chemistry software package. 16

IV. RESULTS AND DISCUSSION
The GUGA-FCIQMC RDM sampling has been used within the Stochastic-CASSCF framework to study the low-energy spin states of the [Fe(III) 2 S 2 (SCH 3 ) 4 ] 2− model complex (Figure 4a), derived from synthetic complexes of Mayerle et al., 123,124 and utilized in our previous investigation, 57 and the [Fe(III) 4 S 4 (SCH 3 ) 44 ] model cubane (Figure 4b), obtained from the synthetic complex of Averill et al., 125 where the terminal groups have been replaced by methyl groups. For the [Fe 2 S 2 ] system, we considered (1) a CAS(10e,10o), consisting of the singly occupied iron 3 d orbitals, (2) a CAS(10e,20o), consisting of the singly occupied iron 3 d and the empty correlating double-shell d orbitals, (3) a CAS(22e,16o) consisting of the singly occupied iron 3 d and the six doubly occupied bridging sulfur 3p orbitals, and (4) a CAS(22e,26o) containing the iron 3 d and d orbitals, and the six bridging sulfur 3 p orbitals. This active space correspond to the one utilized in our previous work. 57 We also studied the role of the iron 4 s and the periph- eral sulfur 3 p orbitals, which were considered in other studies of similar FeS dimers, 95,126,127 having mixedvalence states as main target. We found that the iron 4 s orbitals have a negligible differential role on the lowenergy spin gaps. Including one terminal orbital per peripheral sulfur atom in the active space resulted in an uneven mixing between different orbitals on some of the peripheral sulfur atoms upon completion of the CASSCF procedure. This suggests that for a balanced treatment of the peripheral S orbitals one would need to include all 12 of them. However, while these orbitals have important ligand field effects that could affect the energetic of mixed-valence states, we found that their role is less crucial for dealing with the homo-valent [Fe(III)S] systems. Additionally, a recent study on the excited state spectrum of the [FeS] dimer 128 using the CAS(22e,16o) wave functions, showed that the low-lying non-Hund excited states involve bridging-sulfur charge transfer (CT) states, while CT states involving terminal-sulfur orbitals were only found at higher energies.
Thus, we decided not to further consider these orbitals in the chosen model active space. This was considered to be a successful strategy in previous works. 57,58,129 Similar to previous computational studies 57,58,95,126,127,129 we do not include empty sulfur orbitals in our active space. Therefore, metal-to-ligand charge transfer (MLCT) excitations are not considered by the model active space chosen. However, as also suggested by Neese et al., 127 such configurations are rather high in energy, and they can be safely neglected for low-energy spectrum calculations.
We used an extended relativistic atomic natural orbital basis of double-ζ quality for Fe atoms and a minimal basis for all other elements. The exactly diagonalizable Fe 2 S 2 (10e,10o), (10e,20o) and (22e,16o) active spaces are straightforwardly calculable within spin-adapted Stochastic-CASSCF with modest computational resources. We ensured the convergence of the (22e,26o) active space calculations with respect to the number of walkers, N w , by increasing N w up to N w = 1 · 10 9 , see the SI 97 for more information. The average number of occupied CSFs, N CSF , at each time-step during the GUGA-FCIQMC calculation for N w = 5 · 10 8 and each spin state is shown in Table I

A. Fe2S2 system
The (10e,10o), (10e,20o) and (22e,16o) active spaces are exactly diagonalizable, and were considered to study the differential interplay of different correlation mechanisms, such as orbital relaxation, double-shell 52,130-132 and superexchange 53,133-136 correlation effects, and to benchmark and test our stochastic spin-free RDM sampling procedure. A thorough comparison of the exact and the Stochastic-CASSCF results can be found in the SI 97 .
In our earlier works, 57,58 we have demonstrated via theoretical arguments, and shown with calculations, that the choice of orbital representation and reordering greatly effect the sparsity of the CI wave function within the GUGA formalism. We have also shown that the localization and reordering strategy within the GUGA-FCIQMC algorithm is of utmost importance, as it positively influences the stability of the dynamics and the convergence with respect to the total number of walkers. Moreover, this strategy greatly simplifies the interpretation of the converged wave functions, and could even allow selective optimization of one among groundand low-energy excited-state wave functions. We have adopted the same strategy for the present work. In Reference [57] the optimized CASSCF(22e,26o) orbitals for the S = 0 ground state, obtained via the SD-based Stochastic-CASSCF, 6 were used as starting orbitals for the localization and reordering protocol and for the GUGA-FCIQMC dynamics. A CASSCF(10e,10o) was performed inside the CAS(22e,26o) active space, an invariant rotation within the CAS(22e,26o), that separates valence 3 d orbitals from the six sulfur and the 10 correlating d orbitals. Only the 10 valence 3 d orbitals were localized and site-ordered, leaving the sulfur and the correlating d orbitals delocalized. In the present work, the starting orbitals were obtained from a highspin restricted open-shell Hartree-Fock (ROHF) calcula- tion, equivalent to a CASSCF(10e,10o) S = 5 optimization. The iron 3 d and d orbitals, resulting from the ROHF calculation, were separately localized, using the Pipek-Mezey 137 method, while the bridging sulfur 3 p orbitals were left delocalized. Using localized d orbitals allows to better estimate the local spin of each magnetic center. Fe 2 S 2 spin ladder and total energies Figure 5a shows the spin gaps of all the states relative to the S = 0 ground state -the spin ladder -as function of the total spin after the CASSCF orbital optimization. The spin gaps are lowest in the (10e,10o) active space, with ∆E = 12 mH, between the S = 5 and S = 0 state. The inclusion of the iron d orbital in the (10e,20o) active space qualitatively does not change the obtained spin ladder and it also has a rather smaller quantitative effect, with an only slightly larger ∆E = 17 mH between the S = 5 and S = 0 state. Inclusion of the bridging-sulfur 3 p orbitals has the largest effect on the spin gaps, since it accounts for the metal-bridging ligand correlation, which is differentially more important than the radial correlation effect, 131 accounted for by the inclusion of the d . The consideration of both the iron d and the bridging-sulfur 3 p orbitals in the (22e,26o) active space induces a qualitative change in the obtained spin gaps, which will be further discussed below. Quantitatively the relative spin gaps enlarge by as much as a factor of 3.3, when enlarging the active space, from CAS(10e,10o) to CAS(22e,26o).
In Figure 5b we show the total energy of the S = 0, 3 and 5 states, helpful in describing in absolute terms the correlation effects bound to ligand-to-metal chargetransfer and radial correlation effects. Starting from the CAS(10e,10o), the inclusion of the iron correlating d orbitals, as in the CAS(10e,20o) active space, lowers the total energy more than including the sulfur 3 p orbitals, as in the CAS(22e,16o). The combined inclusion of both iron d and sulfur 3 p orbitals has the surprising effect of lowering the total energies more than the ligand-to-metal charge-transfer and the radial correlation effects on their own. However, the largest differential effect arises from the ligand-to-metal charge-transfer excitations as show in Figure 5a.
Orbital relaxation effect. In this section the overall and the differential effect of the CASSCF orbital relaxation on energies and spin-gaps, together with its effect on the derived model parameters, is discussed. The highest-spin, S = 5, restricted openshell Hartree-Fock (ROHF) orbitals from the (10e,10o) active space are chosen as starting orbitals for all the calculations. The results of the first CASSCF iteration are from here on referred to as CASCI. Figure 6a shows the energy difference of the CASCI results using (10e,10o) ROHF orbitals and the CASSCF results, ∆E = E CASCI − E CASSCF , for the S = 0, 3 and 5 states as a function of the active space. As expected, the effect of the CASSCF orbital relaxation, when using (10e,10o) ROHF orbital, is lowest for the (10e,10o) active (with differences below 10 mH), and highest for the (22e,26o) active space. Within each active space the effect of the CASSCF procedure is largest for the low spin states, with a maximum difference of ∆E = 94 mH for the singlet in the (22e,26o) active space. For the high spin states the effect of the CASSCF procedure is smaller, but still substantial for the larger active spaces, especially in the (22e,26o) AS, with ∆E = 76 mH for the S = 5 state.
To investigate the differential effect, Figure 6b shows the changes in the spin gaps, due to CASSCF orbital relaxation as a function of active space. As expected, the CASSCF procedure increases all the obtained spingaps, as the low spin states are stabilized more by the orbital relaxation when starting from high-spin ROHF orbitals than the higher spin-states, which are better represented by the ROHF orbitals. The change in the spin-gaps is smallest for the (10e,10o), where the ROHF starting orbitals were obtained, and the somehow similar (10e,20o) active space. Interestingly, although the effect of the CASSCF procedure on the total energies is highest for the (22e,26o) active space, see Fig. 6a, the largest effect on the spin gaps is observed in the intermediate (22e,16o) active space. The energy differences to low-spin states, ∆S = 1, 2, 3, are affected only weakly by the CASSCF optimization and stay similar to the CASCI results. This can be explained by the fact that the low-spin states are similarly biased in the ROHF orbital basis, and thus show similar stabilization during the CASSCF procedure. The high-spin states, on the other hand, are less stabilized by the CASSCF procedure, and as a consequence the low-to-high spin-gaps result are enlarged by the orbital relaxation.
There is a substantial differential effect of ≈ 15 − 20 mH due to the CASSCF orbital relaxation, so one has to be cautious when using ROHF orbitals for spin systems, and orbital bias towards the high-spin is to be expected, leading to a systematic under-estimation of spin-gap predictions for anti-ferromagnetically coupled magnetic sites. Even for the seemingly SCF-invariant singlet-triplet spin gap, the CASSCF procedure is crucial to obtain more accurate model magnetic parameters, as it will be discussed below. Figure 7 shows the energy differences with respect to the S = 0 ground state, for the CASCI(22e,26o) (blue circles), and for the CASSCF(22e,26o) (orange square) results. As expected, the spin states are more separated after the CASSCF orbital optimization, with a lowestto-highest spin-state gap nearly doubled by the orbital relaxation effects. Figure 7 also shows the spin ladder ( 6. (a) Change of the total energy for the S = 0, 3 and 5 states and (b) change of the spin gaps relative to the S = 0 state due to the CASSCF orbital relaxation as a function of active space using (10e,10o) ROHF as starting orbitals (CASCI). obtained from mapping the ab initio results to a Heisenberg model, [82][83][84]86,87 with (dashed lines) and without (solid lines) a biquadratic correction. This aspect will be discussed in greater detail in the next section.
Mapping to a spin model As previously done by Sharma et al. 95 and in our laboratories, 58 we map the ab initio low-energy spectrum of the Fe 2 S 2 system to a spin Hamiltonian, as the spinexchange interactions are the dominant form of magnetic interactions in this system. First we map the excitation energies of the Fe 2 S 2 system to the bilinear two-site Heisenberg Hamiltonian with eigenvalues whereŜ A/B are the local spin-5 /2 operators of the two iron centers, and S is the total targeted spin. We obtain the magnetic coupling parameter J by performing a least-squares fit of the energy expression, Equation (16), to the ab initio results of all lowest spin states and study the quality of this mapping as a function of the active space size and the effect of the CASSCF orbital optimization. As shown in Figure 7, the bilinear Heisenberg spin-ladder (solid blue line) models the ab initio CASSCI results with high accuracy. However, minor deviations can be observed for the fitting of the CASSCF results. This finding suggests that orbital relaxation effects account for additional forms of interactions between the metal centers in addition to enlarging the predicted J values. An improved Heisenberg model with biquadratic exchange, [84][85][86][87][88][89][90][91][92] with eigenvalues For the CASSCF results (orange), the extracted J (diamonds) and J (triangles) parameters are larger than the corresponding CASCI results, increasingly so in the larger active spaces, and additionally, the bilinear J and biquadratic J differ. For all but the largest active space, the biquadratic J is about 0.1 mH smaller than the bilinear J, while it is ∼ 0.25 mH larger in the (22e,26o) AS. The differences between the extracted model parameters indicate that a simple bilinear Heisenberg model is not sufficient to describe the CASSCF results.
To quantify this discrepancy and analyze how well a biquadratic model suits the ab initio results we show the relative average error per state ω (in percent) of the corresponding bilinear and biquadratic Heisenberg fits to the CASCI (blue solid and striped) and CASSCF (orange solid and striped) results in Figure 8b . Following Ref. [85], ω is defined as where E C S is the computed ab initio spin gap of spin state S relative to the singlet ground state and E M S is the energy obtained by fitting the bilinear and biquadratic model, Eqs. (16) and (18). N is the number of considered states (with N = 5 in the FeS dimer case, as we only consider the spin-gap relative to the singlet ground state) and ∆E C max is the ab initio energy difference between the S = 5 and singlet state.
The CASCI spin ladders exhibit a clear bilinear Heisenberg behavior, as shown by the small ω values (blue bars) in Fig. 7b. The error is less than 1% for all active space sizes. Larger discrepancies emerge between the CASSCF energies and the bilinear Heisenberg model, indicated by larger ω values (orange striped bars in Fig. 7b). The relative error ω, defined in Eq. (19) takes into account the gap between the S = 0 and S = 5 state, ∆E C max , in the denominator. This causes ω to be largest for the bilinear Heisenberg fit to the CASSCF results in the (10e,10o) active space.
Overall, these discrepancies are still rather small (at most 4%) as shown in Figure 7b, however, not negligible. The biquadratic Heisenberg Hamiltonian, Eq. (17), describes the CASSCF results better, as indicated by a much smaller ω value (less than 1%) in all cases, see Fig. 7b. However, the largest CAS(22e,26o) starts to show deviations also from the biquadratic Heisenberg model. Independently of the quantitative aspects, our calculations confirm the anti-ferromagnetic character of this system, with the CASSCF predicting larger antiferromagnetic magnetic constant than the CASCI procedure. This result, although very promising, is not definitive, and in fact correlation effects, not accounted for in the present work, such as dynamic correlation effects outside the active space and convergence with basis set, could further enhance the deviation from the biquadratic Heisenberg Hamiltonian.
Considering the results of the present work and the ones available in the literature 93,95,96,127,129,138,139 (see the SI 97 for details) some clear trends can be promptly recognized: increasing the active space, performing CASSCF orbital optimization and/or recovering dynamic correlation (MCPDFT), widens the energy spread of the spin ladder, and, thus, yielding a larger effective magnetic coupling coefficient J. The almost doubling of the extracted J and J due to the CASSCF procedure, as seen in Figure 7 and 8a, indicates the important role of orbital relaxation by differentially stabilizing the lowspin state.
This finding clearly shows that one needs to be cautious when using CI energies on ROHF orbitals, and a systematic error is to be expected that overstabilizes higher-spin states over low-spin states. Moreover, the deviation from the simple bilinear Heisenberg model, although small indicates that the complexity of the interactions in [FeS] clusters cannot simply be reduced to a Heisenberg spin-system when aiming at quantitative accuracy; instead, more involved forms of interactions are present, that require complex ab initio Hamiltonians (here exemplified by large CASSCF calculations) and model Hamiltonians (here exemplified by the biquadratic Heisenberg).
CASSCF effects on local spin measurements for Fe 2 S 2 To further investigate the applicability of a (biquadratic) Heisenberg spin model, we look into local spin measurements and spin-spin correlation functions between the two iron centers, and study the CASSCF effect on these quantities. We explain in Appendix A, how we directly measure these quantities and in App. C and D how to extract them from the spin-free 1-and 2-RDMs. We want to emphasize that we are aware that the local spin and spin-spin correlation functions between single and sums of orbitals are representation-dependent quantities. Meaning they are not actual physical observables, but they do depend on the type of employed orbitals, i.e. localized or delocalized orbitals. However, they are still an extremely useful means to provide insight in the chemical and physical properties of compounds and accordingly, are extensively used in the literature. 95,[140][141][142] To ensure reproducibility of our results we want to point out the protocol to obtain the orbitals we used again: The starting orbitals for all calculations, were the (10e,10o) ROHF orbitals, for which the iron 3 d and 3 d were identified and separately localized with the default options of the Pipek-Mezey 137 method in OpenMolcas 15,16 . These orbitals were then relaxed during the Stochastic-CASSCF procedure and the converged orbitals, which remained very localized and in the initial atom-separated order (discussed further below), were used to obtain the corresponding local spin and spin-spin correlation functions. We tested the stability of these results by (a) localizing the final CASSCF orbitals and (b) performing a Procrustes 51,143 transformation to map the CASSCF orbitals as close as possible to the starting ROHF orbitals and found no effect on the obtained local spin and spin-spin correlation functions. Figure 9 shows the local spin expectation value on iron A, Ŝ 2 A , extracted from the CASCI (solid bars), and the CASSCF wave functions (striped bars), for different active spaces and all accessible spin states. The CAS(10e,10o) and CAS(10e,20o) exhibit a local spin expectation value close to the maximum possible, Ŝ 2 A max = 5 2 ( 5 2 + 1) = 8.75, for all spin states. The CASSCF orbital relaxation does not have a significant impact on it. Upon inclusion of the bridging sulfur orbitals in the CAS(22e,16o) -enabling ligandto-metal ("super-exchange-type") excitations -the local spin expectation value remains close to the maximum for CASCI results. However, it substantially drops for all spin states upon CASSCF orbital relaxation. This behavior is enhanced for the CAS(22e,26o), however, for this choice of active space a reduced local spin expectation value for the low-spin states is already obtained for the CASCI calculations. Interestingly, the triplet in the CAS(22e,16o) and CAS(22e,26o), and the quintet in the CAS(22e,26o) have a lower local spin expectation value than the singlet after the CASSCF procedure.
The CASSCF local spin expectation value of the triplet state in the CAS(22e,26o) of Ŝ 2 A min ≈ 6.5 corresponds to a local spin of S A ≈ 2, which raises the question of the applicability of the Heisenberg model mapping. In general for systems with local spin momenta larger than S = 1 /2, local non-Hund excited states, can cause deviations from a pure Heisenberg behavior. 84,87,91,144,145 Additionally, as discussed by Sharma et al. 95 the deviation from the pure S = 5 /2 ion demands accounting for spin and charge delocalization. Both contributions can be related to additional biquadratic terms in the Heisenberg Hamiltonian. [85][86][87]144 One striking advantage of our methodology, based on the FCIQMC algorithm applied onto localized and site-ordered MOs, is that we have direct access to the stochastic representation of the ground state wave function. Thus, to further analyze the deviations from a pure Heisenberg model, we investigated the leading contributions to the CASCI and CASSCF results for each studied active space. As an example, for the CASCI (blue squares) and CASSCF (orange circles) singlet results in the (22e,26o) active space, we show the reference weight (Ref. weight), the sum of all metal-to-metal charge transfer (MMCT), local d → d radial excited configurations (Radial), ligand-to-metal charge transfer (LMCT) and local Hund's rule violating configurations (non-Hund) in a radar plot in Figure 10. It is important to note, that the values are displayed in percent and the radial axes (indicated by the above introduced acronyms) are on a logarithmic scale to allow an easier visual comparison of the different contributions to the ground state wave function and the explicit values can be found in the SI.
The reference weight of the S = 0 state in the (22e,26o) active space drops from a value of 74.4% in the CASCI to 46.1% in the CASSCF wave function. On the other hand, the inter-iron MMCT (F e A 3d ↔ F e B 3d) increase from 6.9% to 12.9% and the bridgingsulfur-to-metal LMCT increase from an already large 13.4% to a substantial 27.9% between the CASCI and CASSCF calculations. Both the radial-type, intra-iron 3d → 3d , (CASCI: 1.5%, CASSCF: 2.1%) and intrairon non-Hund configurations (CASCI: 1.2%, CASSCF: 3.7%) only have marginal contributions in the wave functions. The remaining spin states show similarly large LMCT contributions after the CASSCF procedure in the (22e,26o) active space.
These results suggest that the main driving force in lowering the local spin expectation values are LMCT configurations upon inclusion of the bridging-sulfur orbitals in the active space. However, as shown in Figure 9 by the relatively constant (close-to-maximum) CASCI local spin expectation values for all active space, "just" including the sulfur 3p orbitals does not suffice to correctly capture all relevant correlation mechanisms; instead, the CASSCF orbital relaxation of the ROHF starting orbitals is necessary.
Malrieu et al., 146,147 Angeli and Calzado 148 and Li Manni and Alavi 52 have observed that CASSCF orbitals from a minimal active space are too localized to cor- rectly capture all relevant physical mechanisms in a subsequent second-order multi-reference perturbation theory (MRPT2). This is mainly due to the fact that relevant ligand-to-metal charge transfer (LMCT) excitations do not interact with the zeroth-order wave function due to the generalized Brillouin theorem. [149][150][151] On the other hand, natural magnetic orbitals, obtained by e.g. difference-dedicated CI (DDCI) calculations, [152][153][154][155][156] or optimized CASSCF orbitals from large active space calculations, [52][53][54] show correlation-induced metal-ligand delocalization, by capturing higher-order contributions. 52,146,147 We also studied this effect in the present work by directly comparing the localized high-spin S = 5 (10e,10o) ROHF orbitals (used as the starting orbitals in all CASSCF calculations) with the singlet (22e,26o) CASSCF orbitals. During the Stochastic-CASSCF procedure, performed with OpenMolcas, the orbitals remain quite localized and in the chosen atom-separated order, mentioned above and described in the SI. For reproducibility, it is important to note, that we used the last orbitals of the OpenMolcas CASSCF procedure before the standard final diagonalization of the 1-RDM and transformation to natural (delocalized) orbitals. Furthermore, we performed invariant Procrustes orthogo-nal transformations 51,143 -with the OpenMolcas software package -of the (10e,10o) ROHF iron 3d orbitals to make them as similar as possible to the (22e,26o) singlet CASSCF orbitals, to allow an optimal comparison. Further details of the exact protocol for the comparison and corresponding orbital files can be found in the SI.
In Figure 11 we show the (10e,10o) ROHF (top row) and the CASSCF(22e,26o) singlet (middle row) 3d orbitals of iron A, rendered with the Jmol software package, 157 with an isosurface cutoff value of 0.05. The last row of Figure 11 shows the difference of the corresponding orbitals, computed with the pegamoid.py 158 and Multiwfn software package, 159 and rendered with Jmol with an isosurface cutoff value of 0.007 for all orbitals except the third (3rd column), which has a cutoff value of 0.003 to make differences visible. The delocalization effect of the CASSCF procedure can be seen for orbitals two (2nd column) and four (4th column). The orbital differences show that the CASSCF procedure has a metal-to-ligand delocalization effect, where larger tails of the iron 3d orbitals on the ligands, increases both the kinetic and direct exchange integrals, 134,160 and consequently increasing the absolute value of J. 147,148 As discussed above and shown in Figure 10, the delocalization of the iron 3d orbitals is accompanied by a simultaneous increase of the LMCT contributions in the (22e,26o) CASSCF singlet wave function.
Calzado et al., 147 show a very similar orbital dependence when performing CASCI calculation on extracted J parameters. Their study on local S = 1 binuclear systems shows that the high-spin triplet ROHF orbitals yield a much too low J compared to using singlet or state-specific orbitals. Similarly, Spiller et al., 139 find that when using spin state-averaged CASSCF orbitals, a subsequent NEVPT2 treatment yields lower magnetic coupling than using spin-pure state-specific orbitals. Angeli and Calzado 148 suggest to use average orbitals of the singlet ground and excited states in the minimal active space to include the ionic contributions and thus ligand-metal delocalization and Kubas 128 used spin-averaged Hartree-Fock (SAHF) 161 orbitals for the low-lying excited state spectrum of the [FeS] dimer.
On the other hand, CASSCF misses different physical effects, which tends to emphasize the ionic nature of orbitals 127 and causes MOs of pure ionic wave functions to be too diffuse. 162 Similarly, Malrieu et al. 146 showed that  11. (10e,10o) ROHF (top row) and (22e,26o) S = 0 CASSCF FeA 3d orbitals rendered with Jmol 157 with an isosurface value of 0.05. The difference between the corresponding ROHF and CASSCF orbitals (bottom row) was obtained with Multiwfn 159 and is rendered with an isosurface cutoff of 0.007 (except the 3rd column, which uses a value of 0.003). The protocol to obtain the orbitals and their differences is described in the main text and with more detail in the SI, where also the corresponding orbital files can be found.
the definition of magnetic orbitals from spin-unrestricted DFT calculations strongly overestimate the metal-ligand delocalization, which might be the reason for the rather large J value obtained by BS-DFT 96,138 and DMRG CASCI calculations based on such orbitals. 95 CASSCF effect on spin-spin correlation function for Fe 2 S 2 With a spin-adapted basis and the localized and atom ordered MOs described in the SI, we can use the formulas derived in Appendix B to study the spin-spin interaction between the two magnetic centers in the Fe 2 S 2 system, and the effect of the CASSCF procedure on it. Figure 12 shows the spin-spin correlation function Ŝ A ·Ŝ B between the two magnetic centers from the CASCI (solid), and in the CASSCF wave functions (striped bars), as a function of the active space size for all spin states. For all active spaces the spin-spin alignment changes from anti-ferromagnetic to ferromagnetic, starting from S = 4, as the total spin increases. The spin-spin correlations are somewhat large for the CAS(10e,10o) and CAS(10e,20o), where the CASSCF orbital relaxation does not have a big impact on the expectation values. As for the local spin measurements, the orbital relaxation has the biggest effect in the CAS(22e,26o). The CASSCF procedure has a damping effect on the magnitude of the spin-spin correlations, but does not change the description of the underlying physical behavior of a transition from an anti-ferromagnetic to a ferromagnetic alignment as a function of the total spin. With access to the 1-and 2-RDM we are able to study the spin correlation functions on an orbital-resolved level, including the iron 3d and sulfur 3p orbitals. Figure 13 shows the CASSCF spin-spin correlation function Ŝ 0 ·Ŝ i between the first Fe A 3d orbital and all other orbitals, obtained via the spin-free RDMs. Fig. 13 contains Ŝ 0 ·Ŝ i for all spin states, S = 0 to S = 5 (indicated by the subplot titles) and all active spaces, different color and markers. The x-axes indicate the different orbitals i and different types of orbitals (iron, sulfur, etc.) are separated by vertical dashed lines and data points only show up, when possible. E.g. there are no markers of the (10e,10o) active space results (red triangles) for the iron 4d and sulfur 3p orbitals. The mostly singly oc-cupied first iron A 3d orbital, with index 0, is magnetically parallel aligned to all the other Fe A 3d orbitals, as can be seen by the Ŝ 0 ·Ŝ i ≈ 0.25, ∀i ∈ {Fe A 3d}, all the spin states. Ŝ 0 ·Ŝ i ≈ 1/4 is expected for two ferromagnetically S = 1/2 spins. The magnetic 3d orbitals of iron B are highlighted by the gray background in Fig. 13. Here on can see that with increasing total spin S, indicated by the titles of the subplots, the alignment of the first iron A 3d orbital changes from antiferromagnetic, Ŝ 0 ·Ŝ i < 0, to ferromagnetic alignment, Ŝ 0 ·Ŝ i ≈ 1/4, ∀i ∈ {Fe B 3d}, with the 3d orbitals of iron B. The results confirm that the exchange interaction exclusively happens between the (magnetic) iron 3d orbitals, while the other orbitals are magnetically inert (indicated by a zero value of Ŝ 0 ·Ŝ i ).

B. Fe4S4 system
We now turn to the all-ferric [Fe(III) 4 S 4 (SCH 3 ) 4 )] system. Here we consider the minimal (20e,20o) active space consisting of the iron 3d orbitals of the 4 iron atoms. This active space size is already slightly above the current limit of performing routine FCI calculations. 15,16 Similar to Fe 2 S 2 we performed state-specific and spin-pure Stochastic-CASSCF calculations for all the spin states, from S = 0 up to S = 10. We used the geometry studied in References [58,95], which is, among other computational details, documented in the SI. We used an ANO-RCC-VDZ basis set for Fe and an ANO-RCC-MB for all other elements and ensured that the obtained results are converged w.r.t. number of used walkers N w . They are already with a very modest N w = 1 · 10 6 walkers.
In the Fe 4 S 4 study we use the localized (20e,20o) highspin ROHF orbitals as a starting guess and focus on the effect of the CASSCF orbital relaxation on the extraction of the model parameters and associated physical and chemical interpretations of the results. Similar as in the FeS-dimer section above, we refer to the first iteration of the CASSCF procedure, based on the ROHF orbitals, as CASCI. In our previous work, 58 we found that   the ab initio CASCI results can be very well mapped to a simple bilinear Heisenberg model. However, as seen in Section IV A on the dimer model, orbital relaxation effects can affect the relative energy of the ab initio spin states, and introduce forms of interactions that go beyond the simple bilinear Heisenberg model. Figure 14 shows the energy difference to the S = 0 ground state (markers) and a simple (solid) and a biquadratic (dashed line) Heisenberg fit for CASCI and CASSCF results as a function of the total spin. The CASCI results are well represented by a simple Heisenberg model, whereas the CASSCF results differ from it and necessitate a biquadratic model description, similar to the above studied Fe 2 S 2 case.
To investigate the deviation of the ab initio CASSCF results from a pure Heisenberg model, we computed the local spin and spin-spin correlation for the Fe 4 S 4 system. Figure 15 shows the local spin expectation values of iron A (a), A + B (b) and A + B + C (c) as a function of total spin for the CASCI and CASSCF results. The local spin on the single iron A is close to the maximum possible, (S max A ) 2 = 8.75, for all the spin states, and the effect of the CASSCF orbital relaxation is present, but small. Due to symmetry reasons we can safely assume that this expectation value is equal for all 4 iron centers. The expectation value of the sum of the local spin of the far-distanced irons, A and B, is close to the maximum possible, (S max AB ) 2 = 5(5 + 1) = 30, as can be seen in Figure 15b. This shows that the two far-distanced iron cen- ters are ferromagnetically aligned with S A +S B = 5, as already investigated thoroughly for the singlet ground and excited states in Reference [58]. Again, the orbital relaxation only plays a minor role, but shows the same behavior as for the single iron spin. The local spin expectation value of the sum of three irons, (Ŝ A +Ŝ B +Ŝ C ) 2 , increases from a minimum value close to S min  Figure 15, which is the exact results of a 4-site pure S = 5 /2 Heisenberg model. The close agreement to the theoretically maximal values of the local spin is because we investigate the minimal (20e,20o) active space of the magnetic iron 3d orbitals. Inclusion of ligand orbitals would cause a larger deviation similar to the iron dimer studied above, and as already anticipated in our previous work. 58 We now focus on the spin-spin interaction between the 4 magnetic iron centers. Figure 16 shows the spinspin correlation function, Ŝ i ·Ŝ j between the 4 different iron atoms for the CASCI and CASSCF results as a function of the total spin. Figure 16a confirms that the two iron atoms with the largest distance, A and B, always stay ferromagnetically aligned for all the spin states, with a marginally lowering effect of the CASSCF orbital relaxation. Figure 16b and c show that the spin-spin interaction between two close lying iron atoms, e.g. A − C or A − D, is anti-ferromagnetic for the low spin states, and switches to ferromagnetic alignment for S = 8 and higher. Additionally, these results confirm that these spin-spin interactions are symmetric and that the CASSCF procedure has only a marginal effect on the obtained expectation values.

V. CONCLUSION
In this work, we present our implementation to compute the spin-pure one-and two-body reduced density matrices, via stochastic sampling, within our spin-adapted FCIQMC implementation. This gives us access to spin-pure two-body observables, such as the spin-spin correlation function, and allows to use the GUGA-FCIQMC as a spin-pure CI eigensolver in the spin-pure Stochastic-CASSCF approach (within OpenMolcas). This in turn enables us to stochastically, yet accurately, treat active spaces far larger than conventional CI solvers in a spin-pure manner. The implementation requires only minor modification to the existing GUGA-FCIQMC implementation and introduces only a small computational overhead. This makes the approach quite efficient and allows us to employ up to hundreds of millions of CSFs simultaneously.
We demonstrate the utility of this method by studying two FeS dimer and tetramer model systems. For the dimer, by performing extensive state-specific CASSCF calculations for the lowest state of each accessible spinsymmetry, and four active spaces, we find that: (1) the combined effect of Fe 3d orbital relaxation and the ligand-to-metal charge transfer has a larger influence on the energetics of the spin-ladder than the sum of the two effects alone. (2) When using (10e,10o) ROHF starting orbitals for the CASSCF procedure, its effect is rather small (few mH) on the singlet-triplet gap, while up to ≈ 20 mH for low-spin-high-spin gap. (3) When one maps ab initio results to a (biquadratic) Heisenberg Hamiltonian, performing a spin-pure CASSCF procedure has a large impact on the extracted model parameter.
Access to the spin-pure RDMs with GUGA-FCIQMC, allows us to directly measure local-spin and spin-spin correlation functions. Insight into these quantities, the local (double) occupation number and the electron delocalization effect due to the CASSCF procedure, enable us to argue why the CASCI results using (10e,10o) ROHF orbital agree so well with the bilinear Heisenberg model, while the converged CASSCF do not. The ROHF orbitals are optimized such that the Heisenberg exchange mechanism, is the only possible one. Thus, they are too localized on the iron atoms [146][147][148] and even increasing the active space does not enable to fully capture important spin delocalization and charge fluctuations. 85,95,96 We study the FeS tetramer in the minimal (20e,20o) active space, which in a spin-adapted approach, due to 20 open shell localized 3d orbitals is a formidable task. Also, for the tetramer we find that performing a CASSCF procedure necessitates in the inclusion of the biquadratic term into the spin model to correctly map the ab initio results. GUGA CSFs are not spin-eigenfunctions for intermediate orbitals. It turns out that this ordering, in conjunction with using localized 3d orbitals in the CAS(22e, 26o) for the Fe 2 S 2 system, is even more optimal as the choice studied in Reference 57. More optimal in the sense, that we do have an even higher reference weight (0.67 compared to 0.55 for the singlet CASCI)-more single reference character-and hence a faster convergence. The detailed orbital choice and ordering used can be found in the SI.

Appendix B: Spin-spin interaction
The measurement of local spin quantities additionally allows us to compute the spin-spin interaction between different iron sites.

sites
If we assume a set of local, independent spin operatorŝ S i , with [Ŝ i ,Ŝ j ] = 0 and due to symmetry: Ŝ 2 i = Ŝ 2 j ∀ i, j, we can deduce the following relations: Since, following Eq. (A2), we can measure both (Ŝ A +Ŝ B ) 2 and Ŝ 2 A locally, we can deduce the spin correlation function from purely local spin measurements.
3 sites Next, we consider a model system as the one depicted in Fig. 17 with two long bond distances AB and CD and 4 remaining short distances and assume Ŝ 2 i to be identical for all i. If we assume spin correlations functions to be equal for the same bond distances, e.g. Ŝ A ·Ŝ B = Ŝ C ·Ŝ D and Ŝ A · S C = Ŝ B ·Ŝ C etc, we can deduce, withŜ A +Ŝ B +Ŝ C =Ŝ ABC for short, where again we can measure all quantities on the right of Eq. (B2) directly, via Eq. (A2). S tot for short, we obtain
Appendix C: Orbital resolved local spin and spin correlation function from spin-free RDMs Expressing the local spin operators as 163 with the Pauli matrices 164 σ x = 0 1 1 0 , σ y = 0 −i i 0 , σ z = 1 0 0 −1 (C2) and the fermionic creation (annihilation) operators, a ( †) i,σ of electrons with spin σ in spatial orbital i. This results in the explicit expressions where n iσ = a † iσ a iσ is the fermionic number operator of orbital i and spin σ.
If we express theŜ i ·Ŝ j aŝ and consequently the individual terms aŝ S z i ·Ŝ z j = 1 4 (n i↑ − n i↓ ) (n j↑ − n j↓ ) (C7) S x i ·Ŝ x j = 1 4 a † i↑ a i↓ a † j↑ a j↓ + a † i↑ a i↓ a † j↓ a j↑ + a † i↓ a i↑ a † j↑ a j↓ + a † i↓ a i↑ a † j↓ a j↑ (C8) Combining the x and y terms yieldŝ For i = j we can transform (C10) tô For the total local spin operator this meanŝ and consequently we get the relation With the spin-free excitation operatorsÊ ij = σ a † iσ a iσ , andê ij,kl =Ê ijÊkl − δ jkÊil , we can express (C13) simply as since E ii = σ n iσ andê ii,ii =Ê iiÊii −Ê ii = (n i↑ + n i↓ ) 2 −Ê ii = 2n i↑ n i↓ . With the spin-free excitation operators we can write (C13) as and substituting (C14) on the lhs of (C15), we get leading to the relation allowing us to formulate the apparent spin-dependent quantity S z 2 i entirely in spin-free terms for i = j.
For i = j we can transform (C10) tô With the observation leading to the relation with which the spin-correlation, (C6), can be expressed asŜ i ·Ŝ j = S z i · S z j − 1 2 ê ij,ji + σ n iσ n jσ . (C21) To express (C21) in spin-free terms we can rewrite S z i · S z j − 1 2 σ n iσ n jσ = 1 4 (n i↑ − n i↓ )(n j↑ − n j↓ ) − 1 2 (n i↑ n j↑ + n i↓ n j↓ ) = 1 4 (n i↓ n j↑ − n i↑ n j↓ − n i↓ n j↑ + n i↓ n j↓ ) − 1 2 (n i↑ n j↑ + n i↓ n j↓ ) which allows us to write the spin-spin correlation function entirely in spin-free terms aŝ Eq. (C23) allows us to directly obtain the off-diagonal i = j, spin-spin correlation functions Ŝ i ·Ŝ j from the spin-free 2-RDM elements, ê ij,ji and ê ii,jj on an orbital resolved level.