diff --git a/Makefile b/Makefile
index c7aad79..6cec241 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-IRPF90 = irpf90  --codelet=factor_een_2:100000 #--codelet=factor_een:10000 
+IRPF90 = irpf90  --codelet=factor_een:100000
 FC     = gfortran
 FCFLAGS= -O2 -march=native -ffree-line-length-none -I .
 NINJA  = ninja
diff --git a/el_nuc_el.irp.f b/el_nuc_el.irp.f
index 36ac9ef..038c4fc 100644
--- a/el_nuc_el.irp.f
+++ b/el_nuc_el.irp.f
@@ -1,73 +1,4 @@
-BEGIN_PROVIDER [double precision, factor_een]
- implicit none
- BEGIN_DOC
- ! Electron-electron nucleus contribution to Jastrow factor
- END_DOC
- integer :: i, j, alpha, p, k, l, lmax, cindex
- double precision :: x, y, z, t, c_inv, u, a, b, a2, b2, c, t0
-
- PROVIDE cord_vect
- factor_een = 0.0d0
-
- do alpha = 1, nnuc
-    do j = 1, nelec
-       b = rescale_een_n(j, alpha, 1)
-       do i = 1, nelec
-          u = rescale_een_e(i, j, 1)
-          a = rescale_een_n(i, alpha, 1)
-          a2 = a * a
-          b2 = b * b
-          c = rescale_een_n(i, alpha, 1) * rescale_een_n(j, alpha, 1)
-          c_inv = 1.0d0 / c
-          cindex = 0
-          do p = 2, ncord
-             x = 1.0d0
-             do k = 0, p - 1
-                if ( k /= 0 ) then
-                   lmax = p - k
-                else
-                   lmax = p - k - 2
-                end if
-                t = x
-                do l = 1, rshift(p - k, 1)
-                  t = t * c
-                end do
-                ! We have suppressed this from the following loop:
-                ! if ( iand(p - k - l, 1) == 0 ) then
-                !
-                ! Start from l=0 when p-k is even
-                ! Start from l=1 when p-k is odd
-                if (iand(p - k, 1) == 0) then
-                  y = 1.0d0
-                  z = 1.0d0
-                else
-                  y = a
-                  z = b
-                endif
-                do l = iand(p - k, 1), lmax, 2
-                   ! This can be used in case of a flatten cord_vect
-                   ! cidx = 1 + l + (ncord + 1) * k + (ncord + 1) * (ncord + 1) * (p - 1) + &
-                   !      (ncord + 1) * (ncord + 1) * ncord * (alpha - 1)
-                   cindex = cindex + 1
-                   factor_een = factor_een + cord_vect(cindex, typenuc_arr(alpha)) * (y + z) * t
-                   t = t * c_inv
-                   y = y * a2
-                   z = z * b2
-                end do
-                x = x * u
-             end do
-          end do
-       end do
-    end do
- end do
-
- factor_een = 0.5d0 * factor_een
-
-END_PROVIDER
-
-
-
-BEGIN_PROVIDER [ double precision, factor_een_2 ]
+BEGIN_PROVIDER [ double precision, factor_een ]
  implicit none
  BEGIN_DOC
  !
@@ -76,7 +7,7 @@ BEGIN_PROVIDER [ double precision, factor_een_2 ]
  double precision :: riam, rjam_cn, rial, rjal, rijk
  double precision :: cn
 
- factor_een_2 = 0.0d0
+ factor_een = 0.0d0
 
  do p=2,ncord
    do k=0,p-1
@@ -87,7 +18,9 @@ BEGIN_PROVIDER [ double precision, factor_een_2 ]
      endif
 
      do l=0,lmax
-       if ( iand(p-k-l,1) == 1) cycle
+       if ( iand(p-k-l,1) == 1) then
+         cycle
+       endif
        m = (p-k-l)/2
 
        do a=1, nnuc
@@ -99,7 +32,7 @@ BEGIN_PROVIDER [ double precision, factor_een_2 ]
              rial = rescale_een_n(i,a,l)
              riam = rescale_een_n(i,a,m)
              rijk = rescale_een_e(i,j,k)
-             factor_een_2 = factor_een_2 + &
+             factor_een = factor_een + &
                             rijk * (rial+rjal) * riam * rjam_cn
            enddo
          enddo
@@ -111,3 +44,78 @@ BEGIN_PROVIDER [ double precision, factor_een_2 ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4,nelec) ]
+ implicit none
+ BEGIN_DOC
+ ! dimensions 1-3 : dx,dy,dz
+ !
+ ! Rdimension 4 : d2x + d2y + d2z
+ END_DOC
+ integer :: i,j,a,p,k,l,lmax,m
+ double precision :: riam, rjam_cn, rial, rjal, rijk
+ double precision, dimension(4) :: driam, drjam_cn, drial, drjal, drijk
+ double precision :: cn, v1, v2, l1, l2, d1, d2
+
+ factor_een_deriv_e = 0.0d0
+
+ do p=2,ncord
+   do k=0,p-1
+     if (k /= 0) then
+       lmax = p-k
+     else
+       lmax = p-k-2
+     endif
+
+     do l=0,lmax
+       if ( iand(p-k-l,1) == 1) then
+         cycle
+       endif
+       m = (p-k-l)/2
+
+       do a=1, nnuc
+         cn = cord_vect_lkp(l,k,p,typenuc_arr(a))
+         do j=1, nelec
+           rjal    = rescale_een_n(j,a,l)
+           rjam_cn = rescale_een_n(j,a,m) * cn
+           do ii=1,4
+             drjal(ii)    = rescale_een_n_deriv_e(ii,j,a,l)
+             drjam_cn(ii) = rescale_een_n_deriv_e(ii,j,a,m) * cn
+           enddo
+           factor_een_deriv_e(:,j) = 0.d0
+           do i=1, nelec
+             rial = rescale_een_n(i,a,l)
+             riam = rescale_een_n(i,a,m)
+             rijk = rescale_een_e(i,j,k)
+             do ii=1,4
+               drijk(ii) = rescale_een_e_deriv_e(ii,i,j,k)
+             enddo
+             l1 = 0.d0
+             l2 = 0.d0
+             x(1:3) = 0.d0
+             x(4) = 2.d0
+             do ii=1,4
+               v1 = rijk * (rial+rjal)
+               v2 = rjam_cn * riam
+
+               d1 = drijk(ii) * (rial+rjal) + rijk * (rial+drjal(ii))
+               d2 = drjam_cn(ii) * riam
+
+               l1 = l1 + drijk(ii) * (rial+drjal(ii))
+               l2 = l2 + drjam_cn(ii) * riam
+
+               factor_een_deriv_e(ii,j) = factor_een_deriv_e(ii,j) + &
+                  v1 * d2 + d1 * v2 + x(ii) * (l1 + l2)
+             enddo
+             factor_een_deriv_e(ii,j) = factor_een_deriv_e(ii,j) + &
+                v1 * d2 + d1 * v2
+           enddo
+         enddo
+       enddo
+
+     enddo
+   enddo
+ enddo
+ factor_een_deriv_e = 0.5d0 * factor_een_deriv_e
+
+END_PROVIDER
+
diff --git a/jastrow b/jastrow
index d641b47..a4f3b78 100755
Binary files a/jastrow and b/jastrow differ
diff --git a/jastrow.irp.f b/jastrow.irp.f
index 02d3e71..3b0a3eb 100644
--- a/jastrow.irp.f
+++ b/jastrow.irp.f
@@ -3,6 +3,5 @@ program jastrow
   print *, 'The total Jastrow factor'
   print *, jastrow_full
   print *, factor_een
-  print *, factor_een_2
 
 end program
diff --git a/rescale.irp.f b/rescale.irp.f
index 5fa9c65..2a55d30 100644
--- a/rescale.irp.f
+++ b/rescale.irp.f
@@ -58,10 +58,19 @@ BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec, 0:ncord)]
     enddo
   enddo
  enddo
+ ! More efficient to compute the exp of array than to do it in the loops 
  rescale_een_e = dexp(rescale_een_e)
+
+ ! Later we use a formula looping on i and j=1->j-1. We need to set Rjj=0 to
+ ! enable looping of j=1,nelec do l=0,ncord
+ do l=0,ncord
+   do j=1,nelec
+    rescale_een_e(j, j, l) = 0.d0
+   enddo
+ enddo
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, rescale_een_n, (nelec, nnuc, 0:ncord)]
+BEGIN_PROVIDER [double precision, rescale_een_n, (4, nelec, nnuc, 0:ncord)]
  implicit none
  BEGIN_DOC
  ! R = exp(-kappa r) for electron-electron for $J_{een}$
@@ -79,3 +88,33 @@ BEGIN_PROVIDER [double precision, rescale_een_n, (nelec, nnuc, 0:ncord)]
  enddo
  rescale_een_n = dexp(rescale_een_n)
 END_PROVIDER
+
+BEGIN_PROVIDER [double precision, rescale_een_n_deriv_e, (4,nelec, nnuc, 0:ncord)]
+ implicit none
+ BEGIN_DOC
+ ! R = exp(-kappa r) for electron-electron for $J_{een}$
+ END_DOC
+ integer :: i, j, l
+ double precision :: kappa_l
+
+ do l=0,ncord
+   kappa_l = - dble(l) * kappa
+   do j = 1, nnuc
+     do i = 1, nelec
+       do ii=1,4
+         rescale_een_n_deriv_e(ii, i, j, l) = &
+             kappa_l * elnuc_dist_deriv_e(ii,i,j)
+       enddo
+       rescale_een_n_deriv_e(4, i, j, l) = rescale_een_n_deriv_e(4, i, j, l) + &
+         rescale_een_n_deriv_e(1, i, j, l) * rescale_een_n_deriv_e(1, i, j, l) + &
+         rescale_een_n_deriv_e(2, i, j, l) * rescale_een_n_deriv_e(2, i, j, l) + &
+         rescale_een_n_deriv_e(3, i, j, l) * rescale_een_n_deriv_e(3, i, j, l) 
+       do ii=1,4
+         rescale_een_n_deriv_e(ii, i, j, l) = &
+             rescale_een_n_deriv_e(ii,i,j, l) * rescale_een_n(i, j, l)
+       enddo
+     enddo
+   enddo
+ enddo
+END_PROVIDER
+