pick_random_4ind Subroutine

    subroutine pick_random_4ind(csf_i, elec_1, elec_2, orb_1, orb_2, elecs, orbs, pgen)
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: elec_1, elec_2, orb_1, orb_2
        integer, intent(out) :: elecs(2), orbs(2)
        real(dp), intent(out) :: pgen

        real(dp) :: r

        elecs = 0
        orbs = 0

        pgen = 1.0_dp

        select case (csf_i%stepvector(elec_1))
        case (1)
            elecs(1) = 2 * elec_1 - 1

        case (2)
            elecs(1) = 2 * elec_1

        case (3)
            r = genrand_real2_dSFMT()

            if (r < 0.5_dp) then
                elecs(1) = 2 * elec_1 - 1
            else
                elecs(1) = 2 * elec_1
            end if

            pgen = 0.5_dp

        end select

        select case (csf_i%stepvector(elec_2))
        case (1)
            elecs(2) = 2 * elec_2 - 1

        case (2)
            elecs(2) = 2 * elec_2

        case (3)
            r = genrand_real2_dSFMT()

            if (r < 0.5_dp) then
                elecs(2) = 2 * elec_2 - 1
            else
                elecs(2) = 2 * elec_2
            end if

            pgen = pgen * 0.5_dp

        end select

        select case (csf_i%stepvector(orb_1))
        case (0)
            if (same_spin(elecs(1), elecs(2))) then
                if (is_beta(elecs(1))) then
                    orbs(1) = 2 * orb_1 - 1
                else
                    orbs(1) = 2 * orb_1
                end if
            else
                r = genrand_real2_dSFMT()

                if (r < 0.5_dp) then
                    orbs(1) = 2 * orb_1 - 1
                else
                    orbs(1) = 2 * orb_1
                end if

                pgen = 0.5_dp * pgen
            end if

        case (1)
            orbs(1) = 2 * orb_1

        case (2)
            orbs(1) = 2 * orb_1 - 1

        end select

        ! i think i can restrict the last one or?..

        if (same_spin(elecs(1), elecs(2))) then
            if (is_beta(elecs(1))) then
                if (csf_i%stepvector(orb_2) == 1) then
                    pgen = 0.0_dp
                    elecs = 0
                    orbs = 0
                    return
                else
                    orbs(2) = 2 * orb_2 - 1
                end if
            else
                if (csf_i%stepvector(orb_2) == 2) then
                    pgen = 0.0_dp
                    elecs = 0
                    orbs = 0
                    return
                else
                    orbs(2) = 2 * orb_2
                end if
            end if
        else
            if (is_beta(orbs(1))) then
                if (csf_i%stepvector(orb_2) == 2) then
                    pgen = 0.0_dp
                    elecs = 0
                    orbs = 0
                    return
                else
                    orbs(2) = 2 * orb_2
                end if
            else
                if (csf_i%stepvector(orb_2) == 1) then
                    pgen = 0.0_dp
                    elecs = 0
                    orbs = 0
                    return
                else
                    orbs(2) = 2 * orb_2 - 1
                end if
            end if
        end if

    end subroutine pick_random_4ind