make_triple Subroutine

    subroutine make_triple(nI, nJ, elecs, orbs, ex, tPar)
        integer, intent(in) :: nI(nel), elecs(3), orbs(3)
        integer, intent(out) :: nJ(nel), ex(2, 3)
        logical, intent(out) :: tPar
#ifdef DEBUG_
        character(*), parameter :: this_routine = "make_triple"
#endif
        integer :: sort_elecs(3), sort_orbs(3), src(3), pos_moved, k, i

        ! i should also check if this excitation is actually possible!
        ! and which spin i move to where or??

        ! figure out how to do triples efficiently..
        ! can we do a single and then a double?

        ! NO: talk to Manu and do this with integer representation! not
        ! with nI and nJ, since this can be done way more effective as
        ! via the occupied orbitals..

        ! TODO: thats a super strange convention, .. talk with Ali and
        ! Simon about that.. but for now accept it as it is..

        sort_elecs = sort_unique(elecs)
        sort_orbs = sort_unique(orbs)

        src = nI(sort_elecs)

        ex(1, :) = src
        ex(2, :) = sort_orbs

        nJ = nI

        ASSERT(sum(get_spin_pn(src)) == sum(get_spin_pn(orbs)))

        ! i should do some stuff depending on the order of src and tgt
        ! we have to check if electrons hop over other electrons, so we
        ! might have to change the indexing to adapt to the change in nJ!

        ! 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 (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

        pos_moved = 0

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

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

        end do

        tPar = btest(pos_moved, 0)

    end subroutine make_triple