FindExcitDet Subroutine

    SUBROUTINE FindExcitDet(ExcitMat, nI, IC, TParity)
        integer, intent(in) :: ic
        integer, intent(inout) :: ExcitMat(2, ic)
        integer, intent(inout) :: nI(nel)
        logical, intent(out) :: tParity
#ifdef DEBUG_
        character(*), parameter :: this_routine = "FindExcitDet"
#endif
        integer :: i, k, pos_moved
        integer :: src(ic), sort_orbs(ic), sort_elecs(ic)

    ! adapt this for double excitations now..
#ifdef DEBUG_
        if (IC < 0 .or. ic > 3) then
            call stop_all(this_routine, "wrong ic!")
        end if
#endif

        if (ic == 3 .or. ic == 2) then
    !Ensure i<j and a<b
            ExcitMat(1, :) = sort_unique(ExcitMat(1, :))
            ExcitMat(2, :) = sort_unique(ExcitMat(2, :))
        end if

        sort_elecs = ExcitMat(1, :)

    ! Return Excitmat with orbitals, rather than indices being specified.
        ExcitMat(1, :) = nI(sort_elecs)

    ! this is the same functionality as in make_double.. but much more
    ! involved here.. so just take the core from make double and others..
    ! look there for the code!
        src = nI(sort_elecs)
        sort_orbs = ExcitMat(2, :)

    ! reshuffle the orbitals..
        if (ic > 1) then
    ! or check it individually:
            if (src(2) < sort_orbs(1)) then
    ! then i hops over j:
                sort_elecs(2) = sort_elecs(2) - 1
            end if
            if (ic == 3) then
            if (src(3) < sort_orbs(1)) then
    ! then i hops over k
    ! (note: this also implies that j hops over k, but treat that
    ! seperately below, to also cover the case, where this if here
    ! is not fullfilled!)
                sort_elecs(3) = sort_elecs(3) - 1
            end if
            if (src(3) < sort_orbs(2)) then
    ! then j hops over k
                sort_elecs(3) = sort_elecs(3) - 1
            end if
            end if
        end if

        pos_moved = 0 !This keeps track of how far the common orbitals have moved.

        do k = 1, ic
        if (src(k) < sort_orbs(k)) then
        if (sort_elecs(k) == nel) then
    ! this can only happen for k == 3
            i = nel + 1
            nI(nel) = sort_orbs(k)
        else
            do i = sort_elecs(k) + 1, nel
            if (sort_orbs(k) < nI(i)) then
                nI(i - 1) = sort_orbs(k)
                exit
            else
                nI(i - 1) = nI(i)
            end if
            end do
            if (i == nel + 1) then
                nI(nel) = sort_orbs(k)
            end if
        end if
        else
        if (sort_elecs(k) == 1) then
            i = 0
            nI(1) = sort_orbs(k)
        else
            do i = sort_elecs(k) - 1, 1, -1
            if (sort_orbs(k) > nI(i)) then
                nI(i + 1) = sort_orbs(k)
                exit
            else
                nI(i + 1) = nI(i)
            end if
            end do
            if (i == 0) then
                nI(1) = sort_orbs(k)
            end if
        end if
        end if

        pos_moved = pos_moved + sort_elecs(k) - i + 1

        end do

        tParity = btest(pos_moved, 0)
    END SUBROUTINE FindExcitDet