open_shell_product Function

    function open_shell_product(alpha, beta, n_orbs) result(basis)
        ! compute the tensor product of alpha and beta spin-bases to give
        ! fully open-shell determinants, the first input is to be defined as
        ! the alpha spin part
        ! and the beta input then is only as big as n_orbs - n_alpha and only
        ! tells us in which empty orbitals of the alpha orbitals beta spins
        ! should be set
        ! sorting is not so easily ensured here anymore..
        integer(n_int), intent(in) :: alpha, beta
        integer, intent(in) :: n_orbs
        integer(n_int) :: basis

        integer, allocatable :: nZeros(:), nOnes(:)
        integer(n_int) :: mask_zeros
        integer :: i
        ! the setting of the alpha spins is the same or?
        basis = set_alpha_beta_spins(alpha, n_orbs, .false.)

        ! i need the zeros, but just in the n_orbs range
        mask_zeros = iand(not(alpha), int(maskr(n_orbs), n_int))

        allocate(nZeros(popcnt(mask_zeros)))

        call decode_bit_det(nZeros, [mask_zeros])

        ! now i need the ones from beta
        allocate(nOnes(popcnt(beta)))
        call decode_bit_det(nOnes, [beta])

        ! and i have to set all beta-spins indicated by nZeros(nOnes)
        ! so we need beta spin!
        nZeros = 2 * nZeros - 1
        do i = 1, size(nOnes)
            basis = ibset(basis, nZeros(nOnes(i)) - 1)
        end do

    end function open_shell_product