exchange_contrib_kvec Function

    HElement_t(dp) function exchange_contrib_kvec(nI, J, p, q, k, spin)
        ! the 3-body exchange contribution for the J optimization
        integer, intent(in) :: nI(:), p(N_DIM), q(N_DIM), k(N_DIM), spin
        real(dp), intent(in) :: J
#ifdef DEBUG_
        character(*), parameter :: this_routine = "exchange_contrib_kvec"
#endif
        integer :: k1(N_DIM), k2(N_DIM)

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

        k1 = lat%add_k_vec(p, q)
        k2 = lat%subtract_k_vec(p, q)
        k2 = lat%subtract_k_vec(k2, k)

        exchange_contrib_kvec = -2.0_dp * real(bhub, dp) * (cosh(J) - 1.0_dp) / real(omega, dp) &
                                * (get_one_body_diag(nI, -spin, k1, .true.) + get_one_body_diag(nI, -spin, k2, .true.))

    end function exchange_contrib_kvec