pick_a_orbital_hubbard Subroutine

    subroutine pick_a_orbital_hubbard(ilutI, orb, p_orb, sum_ms)
        ! hm... i think the first orbital can be picked totally random out
        ! of the empty orbitals or? since every spin and every momentum
        ! is allowed.. and i do not want to overdo the weighting i guess,
        ! especially not for the beginning
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: orb
        real(dp), intent(out) :: p_orb
        integer, intent(in), optional :: sum_ms
        integer :: spin

        ! if sum_ms is present, we pick the first orbital from the minority
        ! spins in the picked electrons -> so it is the opposite
        if (present(sum_ms)) then
            if (sum_ms == -1) then
                ! there a are 2 beta and one alpha electron picked ->
                ! so pick alpha here!
                spin = 0
                p_orb = 1.0_dp / real(nbasis / 2 - nOccAlpha, dp)
            else if (sum_ms == 1) then
                spin = -1
                p_orb = 1.0_dp / real(nbasis / 2 - nOccBeta, dp)
            end if

            do
                orb = 2 * (1 + int(genrand_real2_dsfmt() * nbasis / 2)) + spin

                if (IsNotOcc(ilutI, orb)) exit

            end do
        else

            do
                orb = 1 + int(genrand_real2_dsfmt() * nbasis)

                if (IsNotOcc(ilutI, orb)) exit

            end do

            p_orb = 1.0_dp / real(nbasis - nel, dp)
        end if

    end subroutine pick_a_orbital_hubbard