get_one_body_diag_kvec Function

    function get_one_body_diag_kvec(nI, spin, k_shift, t_sign) result(hel)
        integer, intent(in) :: nI(nel)
        integer, intent(in) :: spin, k_shift(3)
        logical, intent(in), optional :: t_sign
        HElement_t(dp) :: hel

        integer :: i, sgn
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_one_body_diag_kvec"
#endif
        ! change this routine to also use just the symmetry symbols
        integer :: sym_shift
        type(symmetry) :: sym
        logical :: t_sign_
        def_default(t_sign_, t_sign, .false.)

        ! the spin input: -1 is beta, +1 is alpha, 0 is both!
        ! if spin is not present, default is both!
        hel = h_cast(0.0_dp)

        ! k_shift is actually always present..
        ! work on the newest, hopefully correct way to do this..
        ! i need -s k vector for the triples contribution to the doubles..
        if (t_sign_) then
            sgn = -1
        else
            sgn = 1
        end if

        sym_shift = lat%k_to_sym(k_shift(1), k_shift(2), k_shift(3))

        if (sgn == 1) then
            ASSERT(spin == -1 .or. spin == 1)
            if (spin == -1) then
                do i = 1, nel
                    if (is_beta(nI(i))) then
                        sym = SymTable(G1(nI(i))%sym%s, sym_shift)
                        hel = hel + epsilon_kvec(sym)
                    end if
                end do
            else if (spin == 1) then
                do i = 1, nel
                    if (is_alpha(nI(i))) then
                        sym = SymTable(G1(nI(i))%sym%s, sym_shift)
                        hel = hel + epsilon_kvec(sym)
                    end if
                end do
            end if
        else if (sgn == -1) then
            ASSERT(spin == -1 .or. spin == 1)
            if (spin == -1) then
                do i = 1, nel
                    if (is_beta(nI(i))) then
                        sym = SymTable(sym_shift, SymConjTab(G1(nI(i))%sym%s))
                        hel = hel + epsilon_kvec(sym)
                    end if
                end do
            else if (spin == 1) then
                do i = 1, nel
                    if (is_alpha(nI(i))) then
                        sym = SymTable(sym_shift, SymConjTab(G1(nI(i))%sym%s))
                        hel = hel + epsilon_kvec(sym)
                    end if
                end do
            end if
        end if

    end function get_one_body_diag_kvec