InitFCIMC_CSF Subroutine

    subroutine InitFCIMC_CSF()
        implicit none
        integer(n_int), allocatable ::  initSpace(:, :)
        integer :: count, nUp, nOpen
        integer :: i, j, lwork, proc
        integer :: DetHash, pos, TotWalkersTmp, nI(nel), nJ(nel)
        integer(n_int) :: ilutJ(0:NIfTot)
        integer(n_int), allocatable :: openSubspace(:)
        real(dp), allocatable :: S2(:, :), eigsImag(:), eigs(:), evs(:, :), void(:, :), work(:)
        real(dp) :: normalization, rawWeight, HDiag, tmpSgn(lenof_sign)
        HElement_t(dp) :: HOffDiag
        integer :: err
        real(dp) :: HFWeight(inum_runs)
        logical :: tSuccess
        character(*), parameter :: t_r = "InitFCIMC_CSF"

        ! get the number of open orbitals
        nOpen = sum(nOccOrbs) - sum(nClosedOrbs)
        ! in a closed shell system, nothing to do
        if (nOpen == 0) then
            call InitFCIMC_HF()
            return
        end if
        ! first, set up the space considered for the CSF
        call generateInitSpace()
        if (allocated(openSubspace)) deallocate(openSubspace)
        ! we now have initSpace(:,:) with iluts belonging to all possible initial
        ! dets (i.e. all dets contributing to the target CSF) -> construct S2
        allocate(S2(count, count))
        do i = 1, count
            do j = 1, count
                S2(i, j) = S2Matel(initSpace(:, i), initSpace(:, j))
            end do
        end do

        ! prepare the diagonalization
        allocate(eigs(count))
        allocate(evs(count, count))
        allocate(eigsImag(count))
        allocate(work(1))
        allocate(void(0, 0))
        ! workspace query, get how much tmp memory we need
        call dgeev('N', 'V', count, S2, count, eigs, eigsImag, void, count, evs, count, work, -1, err)
        ! allocate work array
        lwork = int(work(1))
        deallocate(work)
        allocate(work(lwork))
        ! diagonalize S2
        call dgeev('N', 'V', count, S2, count, eigs, eigsImag, void, count, evs, count, work, lwork, err)
        deallocate(void)
        deallocate(work)
        deallocate(eigsImag)

        ! transfer the eigenvector
        do i = 1, count
            if (abs(S2Init * (S2Init + 1) - eigs(i)) < eps) exit
        end do
        if (i > count) then
            call stop_all(t_r, "Requested S2 eigenvalue does not exist")
        end if
        eigs = evs(:, i)
!      normalization = minval(eigs)
        rawWeight = sum(abs(eigs))
        normalization = InitialPart / rawWeight

!      refLoc = maxloc(eigs)
        eigs = eigs * normalization

        TotWalkers = 0
        iStartFreeSlot = 1
        iEndFreeSlot = 0
        do i = 1, count
            call decode_bit_det(nI, initSpace(:, i))
            proc = DetermineDetNode(nel, nI, 0)
            if (iProcIndex == proc) then
                HDiag = get_diagonal_matel(nI, initSpace(:, i))
                HOffDiag = get_off_diagonal_matel(nI, initSpace(:, i))
                DetHash = FindWalkerHash(nI, size(HashIndex))
                TotWalkersTmp = int(TotWalkers)
                tmpSgn = eigs(i)
                call encode_sign(initSpace(:, i), tmpSgn)
                if (tHPHF) then
                    call FindDetSpinSym(nI, nJ, nel)
                    call encodebitdet(nJ, ilutJ)
                    ! if initSpace(:,i) is not the det of the HPHF pair we are storing,
                    ! skip this - the correct contribution will be stored once
                    ! the spin-flipped version is stored
                    if (DetBitLT(initSpace(:, i), ilutJ, NIfD) == 1) cycle
                end if
                call AddNewHashDet(TotWalkersTmp, initSpace(:, i), DetHash, nI, HDiag, HOffDiag, pos, err)
                TotWalkers = TotWalkersTmp
            end if
            ! reset the reference?
        end do

        call hash_table_lookup(HFDet, ilutHF, nifd, HashIndex, CurrentDets, i, DetHash, tSuccess)
        if (tSuccess) then
            call extract_sign(CurrentDets(:, i), tmpSgn)
            do i = 1, inum_runs
                HFWeight(i) = mag_of_run(tmpSgn, i)
            end do
        else
            HFWeight = 0.0_dp
        end if

        AllTotParts = InitialPart
        AllTotPartsOld = InitialPart
        OldAllAvWalkersCyc = InitialPart
        OldAllHFCyc = HFWeight
        OldAllNoatHF = HFWeight
        if (tNonInitModShift) then
            OldAllNoInitWalk = AllNoInitWalk
            OldAllNoModShiftWalk = AllNoModShiftWalk
            OldAllNoNonModShiftWalk = AllNoNonModShiftWalk
            OldAllNoBosonActiveWalk = AllNoBosonActiveWalk
        end if

        ! cleanup
        deallocate(evs)
        deallocate(eigs)
        deallocate(S2)
        deallocate(initSpace)
    contains

        !------------------------------------------------------------------------------------------!

        subroutine generateOpenOrbIluts()
            use IntegralsData, only: nfrozen
            implicit none

            count = 0
            nUp = (nel + lms) / 2 - sum(nClosedOrbs) + nfrozen / 2
            do i = 1, 2**nOpen - 1
                if (popcnt(i) == nUp) then
                    count = count + 1
                    if (allocated(openSubspace)) openSubspace(count) = i
                end if
            end do
        end subroutine generateOpenOrbIluts

        !------------------------------------------------------------------------------------------!

        subroutine generateInitSpace()
            implicit none

            integer :: nI(nel), nIBase(nel)
            integer :: iEl, iElBase, iOpen
            integer :: openOrbList(nel)
            logical :: previousCont, nextCont, tClosed
            character(*), parameter :: t_r = "generateInitSpace"

            ! create a list of all open-shell determinants with the correct spin+orbs
            nUp = (nel + lms) / 2
            call generateOpenOrbIluts()
            allocate(openSubspace(count))
            call generateOpenOrbIluts()

            ! convert open-shell-only iluts into full iluts
            ! use the reference to determine which orbitals shall participate
            iElBase = 1
            nIBase = 0
            ! generate a list of open orbitals
            iOpen = 0
            openOrbList = 0
            do i = 1, nel
                if (i == 1) then
                    previousCont = .false.
                else
                    previousCont = FDet(i - 1) == FDet(i) - 1
                end if
                if (i == nel) then
                    nextCont = .false.
                else
                    nextCont = FDet(i + 1) == FDet(i) + 1
                end if
                tClosed = .true.
                ! identify open orbitals, using the ordering: does the next one belong
                ! to the same spatial orb?
                if (is_beta(FDet(i)) .and. .not. nextCont) then
                    iOpen = iOpen + 1
                    openOrbList(iOpen) = FDet(i)
                    tClosed = .false.
                end if
                ! or the previous one?
                if (is_alpha(FDet(i)) .and. .not. previousCont) then
                    iOpen = iOpen + 1
                    openOrbList(iOpen) = FDet(i)
                    tClosed = .false.
                end if
                ! if the orbital is not open, it is closed
                if (tClosed) then
                    nIBase(iElBase) = FDet(i)
                    iElBase = iElBase + 1
                end if
            end do
            if (iOpen /= nOpen) then
                write(stderr, *) "nOpen/iOpen conflict", FDet, openOrbList, iOpen, nOpen
                call stop_all(t_r, "Error in determining open shell orbitals")
            end if

            allocate(initSpace(0:NIfTot, count))
            ! now, add the open-shell contribution
            do i = 1, count
                ! start from the closed-shell base
                nI = nIBase
                iEl = iElBase
                ! for each open orb, add the electron
                do j = 1, nOpen
                    ! based on openSubspace(i), we are considering alpha/beta for this
                    ! orb
                    if (btest(openSubspace(i), j - 1)) then
                        nI(iEl) = get_alpha(openOrbList(j))
                    else
                        nI(iEl) = get_beta(openOrbList(j))
                    end if
                    iEl = iEl + 1
                end do
                ! encode the determinant to the initial space
                call sort(nI)
                call EncodeBitDet(nI, initSpace(:, i))
            end do

        end subroutine generateInitSpace

        function S2Matel(ilutA, ilutB) result(matel)
            implicit none
            integer(n_int), intent(in) :: ilutA(0:NIfTot), ilutB(0:NIfTot)
            integer(n_int) :: splus(0:NIfTot), sminus(0:NIfTot)
            real(dp) :: matel
            integer :: k, m, nI(nel), upOrb, downOrb

            matel = 0.0_dp
            if (DetBitEq(ilutA, ilutB, NIfD)) then
                matel = matel + real(lms * (lms + 2), dp) / 4.0_dp
            end if
            ! get the offdiag part of S2: S-S+
            call decode_bit_det(nI, ilutA)
            do k = 1, nel
                ! first, apply S- to all electrons (additively)
                if (is_beta(nI(k))) then
                    sminus = ilutA
                    downOrb = get_alpha(nI(k))
                    clr_orb(sminus, nI(k))
                    set_orb(sminus, downOrb)
                    ! now, apply S+ to all electrons (again, sum)
                    do m = 1, nel
                        ! check if it yields 0
                        if (is_alpha(nI(m)) .or. m == k) then
                            splus = sminus
                            ! if not, go on
                            if (m == k) then
                                clr_orb(splus, downOrb)
                            else
                                clr_orb(splus, nI(m))
                            end if
                            upOrb = get_beta(nI(m))
                            set_orb(splus, upOrb)
                            if (DetBitEq(splus, ilutB, NIFD)) matel = matel + 1.0_dp
                        end if
                    end do
                end if
            end do
        end function S2Matel
    end subroutine InitFCIMC_CSF