three_body_transcorr_fac_kvec Function

    HElement_t(dp) function three_body_transcorr_fac_kvec(nI, p, q, k, spin)
        integer, intent(in) :: nI(nel), p(N_DIM), q(N_DIM), k(N_DIM), spin
#ifdef DEBUG_
        character(*), parameter :: this_routine = "three_body_transcorr_fac_kvec"
#endif
        integer :: k1(3), k2(3)

        ASSERT(spin == 1 .or. spin == -1)

        ! update: add k-vec to not leave first BZ
        k1 = lat%subtract_k_vec(k, q)
        k2 = lat%add_k_vec(p, q)
        three_body_transcorr_fac_kvec = -three_body_prefac * ( &
                                        n_opp(spin) * (epsilon_kvec(p) + epsilon_kvec(k)) - ( &
                                        get_one_body_diag(nI, -spin, k1) + get_one_body_diag(nI, -spin, k2, .true.)))

    end function three_body_transcorr_fac_kvec