get_diag_helement_k_sp_hub Function

    function get_diag_helement_k_sp_hub(nI) result(hel)
        integer, intent(in) :: nI(nel)
        HElement_t(dp) :: hel
        integer :: i, j, id(nel), idX, idN, k
        HElement_t(dp) :: hel_sing, hel_doub, hel_one, hel_three
        HElement_t(dp) :: temp_hel
        type(symmetry) :: p_sym, k_sym

        ! todo: in the case of 2-body-transcorrelation, there are more
        ! contributions..
        hel = h_cast(0.0_dp)

        if (t_trans_corr_2body) then
            ! this is the
            hel_sing = sum(GetTMatEl(nI, nI))

            id = get_spatial(nI)

            hel_doub = h_cast(0.0_dp)
            hel_one = h_cast(0.0_dp)
            hel_three = h_cast(0.0_dp)

            temp_hel = 0.0_dp

            ! redo this whole shabang.. the formulas are actually much easier:
            ! but just to be sure for now, do i explicetly without any use of
            ! symmetry
            do i = 1, nel
                do j = 1, nel
                    ! the restriction is, that i and j must have opposite
                    ! spin! this also excludes i == j
                    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
                        hel_doub = hel_doub + 0.5_dp * get_umat_kspace(idN, idX, idN, idX)

                        ! then we need the factor of the one-body transcorr influence
                        ! t is defined as -t in our code!, so bhub is usually -1
                        ! and look in the formulas it is actually -2t*cos(k)*2(cosh J - 1)
                        ! (with the k-vector of orbial i!
                        hel_one = hel_one + epsilon_kvec(G1(nI(i))%Sym) &
                                  * omega * three_body_prefac

                        ! and the next part is the three-body with the direct
                        ! and the exchange parts
                        do k = 1, nel
                            ! and the convention is that j and k have same spin!
                            ! and j == k is also allowed and part of it!
                            if (j == k) cycle
                            if (same_spin(ni(j), nI(k))) then
                                ! the k vector is of i and i + j - k
                                ! i need the electrons here ofc..
                                ! even better then the correct k-vector addition
                                ! would be to store an epsilon-k in terms of
                                ! the symmetry symbols!
                                ! something like this but nicer!
                                p_sym = G1(nI(i))%sym
                                k_sym = SymTable(G1(nI(j))%sym%s, SymConjTab(G1(nI(k))%sym%s))

                                hel_three = hel_three - three_body_prefac * ( &
                                            epsilon_kvec(p_sym) - &
                                            (epsilon_kvec(SymTable(p_sym%s, k_sym%s))))

                            end if
                        end do
                    end if
                end do
            end do

            hel = hel_sing + hel_doub + hel_one + hel_three

        else
            hel = sltcnd_excit(nI, Excite_0_t())
        end if

    end function get_diag_helement_k_sp_hub