create_heisenberg_fock_space Function

    function create_heisenberg_fock_space(n_orbs) result(heisenberg_fock)
        ! this routine creates the full heisenberg Fock space for a
        ! given number of orbitals. this can be used in the
        ! calculation of orbital reduced density matrices
        integer, intent(in) :: n_orbs
        integer(n_int), allocatable :: heisenberg_fock(:,:)

        integer :: i, n_beta, n_alpha, j
        integer(n_int), allocatable :: hilbert_space(:,:)

        allocate(heisenberg_fock(0:0, 2**n_orbs), source = 0_n_int)

        j = 1
        do i = 0, n_orbs
            n_beta = n_orbs - i
            n_alpha = i

            hilbert_space = create_all_open_shell_dets(n_orbs, n_beta, n_alpha)
            heisenberg_fock(:,j:j+size(hilbert_space,2)-1) = hilbert_space

            j = j + size(hilbert_space,2)
        end do

    end function create_heisenberg_fock_space