# get_diag_helement_k_sp_hub Function

## public function get_diag_helement_k_sp_hub(nI) result(hel)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)

## Source Code

    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