FindBitExcitLevel_hphf Function

    pure function FindBitExcitLevel_hphf(ilutnI, iLutnJ) result(ic)
        integer(n_int), intent(in) :: ilutnI(0:nifd), ilutnJ(0:nifd)
        integer :: ic
        integer(n_int) :: ilutsI(0:nifd, 2), ilutsJ(0:nifd, 2), tmp(0:nifd)
        integer :: ic_tmp(4), i, j, k

        ilutsI(:, 1) = ilutnI
        call spin_sym_ilut(ilutnI, ilutsI(:, 2))
!         ilutsI(:,2) = spin_flip(ilutnI)

        ilutsJ(:, 1) = ilutnJ
        call spin_sym_ilut(ilutnJ, ilutsJ(:, 2))
!         ilutsJ(:,2) = spin_flip(ilutnJ)

        i = 1
        ic_tmp = 9999

        do j = 1, 2
            do k = 1, 2
                tmp = ieor(ilutsI(:, j), ilutsJ(:, k))
                tmp = iand(ilutsI(:, j), tmp)

                ic_tmp(i) = CountBits(tmp, nifd)
                i = i + 1
            end do
        end do

        ic = minval(ic_tmp)

    end function FindBitExcitLevel_hphf