get_2_body_diag_transcorr Function

    function get_2_body_diag_transcorr(nI) result(two_body)
        integer, intent(in) :: nI(nel)
        HElement_t(dp) :: two_body
        integer :: i, j, id(nel), idX, idN

        two_body = h_cast(0.0_dp)

        id = get_spatial(nI)

        do i = 1, nel
            do j = 1, nel
                if (.not. same_spin(nI(i), nI(j))) then

                    idX = max(id(i), id(j))
                    idN = min(id(i), id(j))

                    ! now we need 1/2, since we loop over all electrons
                    two_body = two_body + 0.5_dp * get_umat_kspace(idN, idX, idN, idX)

                    two_body = two_body + epsilon_kvec(G1(nI(i))%Sym) &
                              * omega * three_body_prefac

                end if
            end do
        end do

    end function get_2_body_diag_transcorr