contract_2_rdm_ind Function

    pure function contract_2_rdm_ind(i, j, k, l, excit_lvl, excit_typ) result(ijkl)
        ! since I only ever have spatial orbitals in the GUGA-RDM make
        ! the definition of the RDM-index combination differently
        integer, intent(in) :: i, j, k, l
        integer, intent(in), optional :: excit_lvl, excit_typ
        integer(int_rdm) :: ijkl

        integer(int_rdm) :: ij, kl

        ij = contract_1_rdm_ind(i, j)
        kl = contract_1_rdm_ind(k, l)

        ijkl = (ij - 1) * (nSpatOrbs**2) + kl

        if (present(excit_lvl)) then
            call encode_excit_lvl_rdm(ijkl, excit_lvl)
        end if

        if (present(excit_typ)) then
            call encode_excit_typ_rdm(ijkl, excit_typ)
        end if

    end function contract_2_rdm_ind