get_one_body_diag_sym Function

    function get_one_body_diag_sym(nI, spin, k_sym, t_sign) result(hel)
        integer, intent(in) :: nI(nel)
        integer, intent(in) :: spin
        type(symmetry), intent(in) :: k_sym
        logical, intent(in), optional :: t_sign
        HElement_t(dp) :: hel

        integer :: i, sgn
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_one_body_diag_sym"
#endif
        ! change this routine to also use just the symmetry symbols
        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_sym 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

        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, k_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(G1(nI(i))%sym%s, k_sym%s)
                        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(k_sym%s, 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(k_sym%s, 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_sym