  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.
