get_ind Function

    function get_ind(this, elec_ind, orb_ind) result(ind)
        ! depending on the cc_order this encodes all the electron and
        ! orbital indices of an excitation in a unique linear fashion
        ! we can assume the electron and orbital indices to be ordered from
        ! highest to lowest
        class(cc_amplitude), intent(in) :: this
        integer, intent(in) :: elec_ind(this%order), orb_ind(this%order)
        integer :: ind

        character(*), parameter :: this_routine = "get_ind"

        integer :: ij, ab, nij, orb(2)

        ind = -1

        ! to i want to access this with invalid excitations?
        ASSERT(all(elec_ind > 0))
        ASSERT(all(orb_ind <= nbasis))

        ! and assert ordered inputs..
        ! or do we want to order it here, if it is not ordered?
        ! tbd!
        ASSERT(minloc(elec_ind, 1) == 1)
        ASSERT(maxloc(elec_ind, 1) == this%order)

        ASSERT(minloc(orb_ind, 1) == 1)
        ASSERT(maxloc(orb_ind, 1) == this%order)

        select case (this%order)
        case (1)
            ! single excitation
            ! the elec_ind goes from 1 to nel and the orb_ind goes from
            ! nel + 1 to nbasis (has nbasis - nel) possible values
            ! can we assume a closed shell ordered reference?
            ! with the nel lowest orbitals occupied?
            ! otherwise this is a bit more tricky..
            ind = ind_matrix_singles(elec_ind(1), orb_ind(1))
            return

            ind = (elec_ind(1) - 1) * (nbasis - nel) + (orb_ind(1) - nel)

        case (2)

            ij = linear_elec_ind(elec_ind(1), elec_ind(2))
            ab = linear_orb_ind(orb_ind(1), orb_ind(2))

            if (ij == 0 .or. ab == 0) then
                ind = 0
            else
                ind = ind_matrix_doubles(ij, ab)
            end if

            return

            call stop_all(this_routine, "still buggy for double excitations!")
            ! double excitation
            ! first encode the ij electron index in a linear fashion
            ASSERT(elec_ind(1) < elec_ind(2))
            ASSERT(orb_ind(1) < orb_ind(2))

            orb = orb_ind - nel

            ij = (elec_ind(1) - 1) * nel - elec_ind(1) * (elec_ind(1) - 1) / 2 + (elec_ind(2) - elec_ind(1))
!             ij = (elec_ind(2) - 1)*elec_ind(2)/2 + elec_ind(1) - 1
            ! i shift the orb_indices by nel..
!             ab = (orb_ind(1) - nel - 1)*(nbasis - nel) + (orb_ind(2) - orb_ind(1))
            ab = (orb(1) - 1) * (nbasis - nel) - orb(1) * (orb(1) - 1) / 2 + (orb(2) - orb(1))
!             ab = (orb_ind(2) - nel - 1)*(orb_ind(2)-nel)/2 + orb_ind(1) - nel -1

            nij = nel * (nel - 1) / 2
            ind = (ij - 1) * nij + ab

        case default
            call stop_all(this_routine, "higher than cc_order 2 not yet implemented!")
        end select

    end function get_ind