get_ex Function

    function get_ex(this, ind) result(ex)
        ! this function gives the specific excitation, if present in the
        ! cc_amplitudes, which are encoded in a linear fashion
        ! with the convention that all the electron and orbital indices are
        ! always provided in an ordered(from lowest to highest) fashion
        class(cc_amplitude), intent(in) :: this
        integer, intent(in) :: ind
        integer :: ex(2, this%order)
        character(*), parameter :: this_routine = "get_ex"

        integer :: ij, ab, cum, i, j, a, b, nij

        ASSERT(ind > 0)
        select case (this%order)
        case (1)

            if (allocated(this%operators)) then
                ex = this%operators(ind, :, :)
                return
            end if

            ! fix this below if i want to actually encode that directly and
            ! not store the operators..
            ASSERT(ind <= nel * (nbasis - nel))

            ! it is a single excition so this is easy to decompose
            ! this decomposition depends on the encoding below!
            ex(2, 1) = mod(ind - 1, (nbasis - nel))
            ! lets hope the integer division is done correctly.. on all compilers
            ex(1, 1) = int((ind - ex(2, 1)) / (nbasis - nel)) + 1
            ! and modify by nel the orbital index again to get the real
            ! orbital index!
            ex(2, 1) = ex(2, 1) + nel + 1

        case (2)

            if (allocated(this%operators)) then
                ex = this%operators(ind, :, :)
                return
            end if

            ! fix this below if i want to actually encode that directly and
            ! not store the operators..
            call stop_all(this_routine, "still buggy for double excitations!")
            ! first we have to get ij, ab back:

            nij = nel * (nel - 1) / 2

            ab = mod(ind - 1, nij)
            ij = (ind - ab) / (nij) + 1

            ab = ab + 1

            ! then we have to get ij and ab from those indices in the same
            ! way
            ! but here it gets more tricky because we can just do the mod..
            ! maybe i have to store the indices in the end..
            i = 1
            cum = (nel - i)

            ! do a stupid sum..
            do while (cum < ij .and. i <= nel)
                i = i + 1
                cum = cum + (nel - i)
            end do

            j = ij + i - (cum - nel + i)

            ! and the same for ab
            a = 1
            cum = (nbasis - nel - a)
            do while (cum < ab .and. a <= (nbasis - nel))
                a = a + 1
                cum = cum + (nbasis - nel - a)
            end do

            b = ab + a - (cum - (nbasis - nel) + a)

            ex(1, :) = [i, j]
            ex(2, :) = [a, b] + nel

        case (3)
            ex = this%operators(ind, :, :)

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

        end select

    end function get_ex