rpa_contrib_kvec Function

    HElement_t(dp) function rpa_contrib_kvec(J, p, k, spin)
        ! gives the rpa contribution in the J-optimization
        real(dp), intent(in) :: J
        integer, intent(in) :: p(N_DIM), k(N_DIM), spin
#ifdef DEBUG_
        character(*), parameter :: this_routine = "rpa_contrib_kvec"
#endif
        integer :: q(N_DIM)

        q = lat%subtract_k_vec(p, k)

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

        rpa_contrib_kvec = real(bhub, dp) * (cosh(J) - 1.0_dp) / real(omega, dp) * &
                           (n_opp(spin) - 1.0_dp) * (epsilon_kvec(p) + epsilon_kvec(q))

    end function rpa_contrib_kvec