create_all_spin_flips Function

    function create_all_spin_flips(nI_in) result(spin_flips)
        ! takes a given spin-configuration in nI representation and
        ! creates all possible states with flipped spin and same ms
        integer, intent(in) :: nI_in(:)
        integer, allocatable :: spin_flips(:, :)

        integer :: nI(size(nI_in)), nJ(size(nI_in))
        integer :: num, n_open, ms, n_states, i, j, k, n_found
        integer(n_int) :: ilutI(0:niftot), ilutJ(0:niftot)
        integer, allocatable :: open_shells(:)
        integer(n_int), allocatable :: ilut_list(:, :)

        nI = nI_in

        num = size(nI)

        ! it is not necessary to have nI sorted at input, so do that now
        call sort_orbitals(nI)

        call EncodeBitDet(nI, ilutI)

        n_open = count_open_orbs(ilutI)
        ms = sum(get_spin_pn(nI))

        ! the number of possible spin distributions:
        n_states = int(choose_i64(n_open, n_open / 2 + ms))

        allocate(spin_flips(num, n_states))
        spin_flips = 0
        ! the first det will be the original one
        spin_flips(:, 1) = nI

        allocate(ilut_list(0:niftot, n_states))
        ilut_list = 0_n_int
        ilut_list(:, 1) = ilutI

        n_found = 1
        ! i have a brute force solution for now:
        ! loop over all the possible states, starting with the original,
        ! since we need to account for multiple flips in this way
        ! i dont need to take the last one..
        do i = 1, n_states - 1
            ! here determine the positions of the open-shell orbitals
            open_shells = find_open_shell_indices(spin_flips(:, i))

            ! and now flip all possible open-shell pairs
            do j = 1, n_open - 1
                do k = j + 1, n_open
                    ! cycle if same orbital or parallel spin
!                     if (j == k) cycle
                    if (same_spin(spin_flips(open_shells(j), i), spin_flips(open_shells(k), i))) cycle

                    nj = spin_flips(:, i)
                    if (is_beta(spin_flips(open_shells(j), i))) then
                        ! then we know (j) is beta and (k) is alpha
                        nJ(open_shells(j)) = &
                            get_alpha(spin_flips(open_shells(j), i))
                        nJ(open_shells(k)) = &
                            get_beta(spin_flips(open_shells(k), i))
                    else
                        ! otherwise (j) is alpha and (k) is beta
                        nJ(open_shells(j)) = &
                            get_beta(spin_flips(open_shells(j), i))
                        nJ(open_shells(k)) = &
                            get_alpha(spin_flips(open_shells(k), i))
                    end if

                    ! if this determinant is not yet in the spin-flip
                    ! list, put it in
                    call EncodeBitDet(nJ, ilutJ)
                    if (.not. is_in_list_ilut(ilutJ, n_found, ilut_list)) then
                        n_found = n_found + 1
                        ilut_list(:, n_found) = ilutJ
                        spin_flips(:, n_found) = nJ
                    end if
                end do
            end do
        end do

!         call print_matrix(transpose(spin_flips))

    end function create_all_spin_flips