FindBitExcitLevel Function

    pure function FindBitExcitLevel(iLutnI, iLutnJ, maxExLevel, t_hphf_ic) result(IC)

        ! Find the excitation level of one determinant relative to another
        ! given their bit strings (the number of orbitals they differ by)
        !
        ! In:  iLutnI, iLutnJ    - The bit representations
        !      maxExLevel        - An (optional) maximum ex level to consider
        !      t_hphf_ic         - An (optional) flag to determine the
        !                          minimum excitation level in an HPHF calculation to
        !                          both spin-coupled references if present
        ! Ret: FindBitExcitLevel - The number of orbitals i,j differ by

        integer(kind=n_int), intent(in) :: iLutnI(0:NIfD), iLutnJ(0:NIfD)
        integer, intent(in), optional :: maxExLevel
        logical, intent(in), optional :: t_hphf_ic
        integer(kind=n_int) :: tmp(0:NIfD)
        integer :: IC, unused

        ! Unused
        if (present(maxExLevel)) unused = maxExLevel

        if (present(t_hphf_ic)) then
            if (t_hphf_ic .and. tHPHF) then
                if (.not. (TestClosedShellDet(ilutnI) .and. TestClosedShellDet(iLutnJ))) then
                    ! make sure that we are calculating the correct excitation
                    ! level, which should be the minimum of the possible ones in
                    ! HPHF mode
                    ! if both are closed shell it is fine
                    ic = FindBitExcitLevel_hphf(ilutnI, ilutnJ)
                    return
                end if
            end if
        end if

        ! Obtain a bit string with only the excited orbitals one one det.
        tmp = ieor(iLutnI, iLutnJ)
        tmp = iand(iLutnI, tmp)

        IC = CountBits(tmp, NIfD)

    end function FindBitExcitLevel