conjugate_rdm_ind_vec Function

    pure function conjugate_rdm_ind_vec(rdm_ind, order) result(conj_rdm_ind)
        integer(int_rdm), intent(in) :: rdm_ind(:)
        integer, intent(in) :: order
        integer(int_rdm), allocatable :: conj_rdm_ind(:)

        integer :: i, j, k, l, m

        allocate(conj_rdm_ind(size(rdm_ind)), source=0_int_rdm)

        if (order == 1) then

            do m = 1, size(rdm_ind)
                call extract_1_rdm_ind(rdm_ind(m), i, j)
                conj_rdm_ind(m) = contract_1_rdm_ind(j, i)
            end do
        else if (order == 2) then
            do m = 1, size(rdm_ind)
                call extract_2_rdm_ind(rdm_ind(m), i, j, k, l)
                conj_rdm_ind(m) = contract_2_rdm_ind(k, l, i, j)
            end do
        end if

    end function conjugate_rdm_ind_vec