A new method to perform complete active space second-order perturbation theory on top of large active spaces optimized with full configuration quantum Monte Carlo is presented. Computing the three- and Fock-contracted four-particle density matrix from imaginary-time-averaged wave functions is found to resolve fermionic positivity violations and to ensure numerical stability. The protocol is applied to [NiFe]-hydrogenase, [Cu2O2]-oxidase and Fe-porphyrin model systems up to 26 electrons in 27 orbitals and benchmarked against DMRG-CASPT2.