setup_ind_matrix_doubles Subroutine

    subroutine setup_ind_matrix_doubles()
        character(*), parameter :: this_routine = "setup_ind_matrix_doubles"

        integer :: n_elec_pairs, i, j, a, b, elec_i, elec_j
        integer :: ij, ab, orb_a, orb_b, k
        logical :: t_par
        integer(n_int) :: temp_ilut(0:nifd)

        ASSERT(allocated(projedet))
        ASSERT(allocated(iLutRef))

        nvirt = nbasis - nel
        n_virt_pairs = nvirt * (nvirt - 1) / 2
        n_elec_pairs = nel * (nel - 1) / 2

        call setup_virtual_mask()

        call setup_elec_ind_mat()
        call setup_orb_ind_mat()

        ! encode only the possible excitations from the closed shell
        ! reference!
        ! i < j; a < b and (i,j) /= (a,b)

        ASSERT(n_elec_pairs == ElecPairs)

        allocate(ind_matrix_doubles(n_elec_pairs, n_virt_pairs))
        ind_matrix_doubles = 0

        temp_ilut = iLutRef(0:nifd, 1)
        k = 1

        ! i could also just loop over the reference determinant..
        do i = 1, nel
            elec_i = projedet(i, 1)
            do j = i + 1, nel
                elec_j = projedet(j, 1)

                t_par = same_spin(elec_i, elec_j)

                ! i have to convert (i,j) -> to a linear index
                ! and i definetly have to use the actual spin orbitals,
                ! since this is the quantity we have access to in the rest of
                ! the calculation
                ij = linear_elec_ind(elec_i, elec_j)
                ! maybe here i could to a check if ij = 0 ?

                ASSERT(ij > 0)
                ASSERT(ij <= ElecPairs)

                ! now i have to check if the orbitals fit and if the spin
                ! fits!
                ! use the virtual mask from back-spawn!
                if (t_par) then
                    do a = 1, Nvirt
                        orb_a = mask_virt_ni(a, 1)
                        do b = a + 1, Nvirt
                            orb_b = mask_virt_ni(b, 1)

                            if (same_spin(orb_a, orb_b) .and. same_spin(orb_a, elec_i)) then
                                ab = linear_orb_ind(orb_a, orb_b)
                                ASSERT(ab > 0)
                                ASSERT(ab <= n_virt_pairs)

                                ind_matrix_doubles(ij, ab) = k
                                k = k + 1

                            end if
                        end do
                    end do
                else
                    do a = 1, nvirt
                        orb_a = mask_virt_ni(a, 1)
                        do b = a + 1, nvirt
                            orb_b = mask_virt_ni(b, 1)
                            if (.not. same_spin(orb_a, orb_b)) then
                                ab = linear_orb_ind(orb_a, orb_b)

                                ASSERT(ab > 0)
                                ASSERT(ab <= n_virt_pairs)

                                ind_matrix_doubles(ij, ab) = k
                                k = k + 1
                            end if
                        end do
                    end do
                end if
            end do
        end do

    end subroutine setup_ind_matrix_doubles