! ---



subroutine jBH_elem_fct_grad(alpha, r1, r2, fct, grad1_fct)

  implicit none
  double precision, intent(in)  :: alpha, r1(3), r2(3)
  double precision, intent(out) :: fct, grad1_fct(3)
  double precision              :: dist, tmp1, tmp2, dist_inv

  dist =   (r1(1) - r2(1)) * (r1(1) - r2(1)) &
         + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
         + (r1(3) - r2(3)) * (r1(3) - r2(3))


  if(dist .ge. 1d-15) then
    dist_inv = 1.d0/dsqrt( dist )
    dist = dist_inv * dist

    tmp1 = 1.d0 / (1.d0 + alpha * dist)

    fct = alpha * dist * tmp1
    tmp2 = alpha * tmp1 * tmp1 * dist_inv
    grad1_fct(1) = tmp2 * (r1(1) - r2(1))
    grad1_fct(2) = tmp2 * (r1(2) - r2(2))
    grad1_fct(3) = tmp2 * (r1(3) - r2(3))

  else

    grad1_fct(1) = 0.d0
    grad1_fct(2) = 0.d0
    grad1_fct(3) = 0.d0
    fct = 0.d0

  endif

  return
end

! ---