get_3_body_diag_transcorr Function

    function get_3_body_diag_transcorr(nI) result(three_body)
        integer, intent(in) :: nI(nel)
        HElement_t(dp) :: three_body

        integer :: i, j, k
        type(symmetry) :: p_sym, k_sym

        three_body = h_cast(0.0_dp)

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

                    do k = 1, nel

                        if (j == k) cycle

                        if (same_spin(nI(j), nI(k))) then
                            p_sym = G1(nI(i))%sym
                            k_sym = SymTable(G1(nI(j))%sym%s, SymConjTab(G1(nI(k))%sym%s))

                            three_body = three_body - 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

    end function get_3_body_diag_transcorr