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.

19 pages, 17 Figures