diff --git a/Makefile b/Makefile
index c94629d..5de1932 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,8 @@
-IRPF90 = irpf90 -a -d
-FC     = gfortran
-FCFLAGS= -O2 -ffree-line-length-none -I .
+IRPF90 = irpf90 --codelet factor_een:100
+#FC = gfortran
+#FCFLAGS= -g -msse4.2  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
+FC     = ifort -g
+FCFLAGS= -O2 -xHost -I .
 NINJA  = ninja
 AR     = ar
 RANLIB = ranlib
diff --git a/el_nuc_el.irp.f b/el_nuc_el.irp.f
index 081693a..ae735e5 100644
--- a/el_nuc_el.irp.f
+++ b/el_nuc_el.irp.f
@@ -4,24 +4,33 @@ BEGIN_PROVIDER [double precision, factor_een]
  ! Electron-electron nucleus contribution to Jastrow factor
  END_DOC
  integer :: i, j, alpha, p, k, l, lmax = 0
+ double precision :: x, y, z, t, t_inv
  factor_een = 0.0d0
  
  do alpha = 1, nnuc
     do j = 1, nelec
        do i = 1, nelec
+        t_inv = 1.d0/(rescale_een_n(i, alpha) * rescale_een_n(j, alpha))
         do p = 2, ncord
-           do k = p - 1, 0
+           do k = p - 1, 0, -1
               if ( k == 0 ) then
                  lmax = p - k - 2
               else
                  lmax = p - k
               end if
-              do l = lmax, 0
-                 if ( mod(p - k - l, 2) == 0 ) then
-                    factor_een = factor_een + cord_vect(p, k, l) * rescale_een_e(i, j) ** k &
-                         * (rescale_een_n(i, alpha) ** l + rescale_een_n(j, alpha) ** l) *  &
-                         (rescale_een_n(i, alpha) * rescale_een_n(j, alpha)) ** ((p - k - l) * 0.5d0)
+              x = 1.d0
+              y = 1.d0
+              z = 1.d0
+              t = (rescale_een_n(i, alpha) * rescale_een_n(j, alpha)) ** (rshift(p - k,1))
+              do l = 0, lmax
+                 if ( iand(p - k - l, 1) == 0 ) then
+                    factor_een = factor_een + cord_vect(l, k, p, alpha) * x &
+                         * (y + z) * t
+                    t = t * t_inv
                  end if
+                 x = x * rescale_een_e(i, j)
+                 y = y * rescale_een_n(i, alpha) 
+                 z = z * rescale_een_n(j, alpha) 
               end do
            end do
         end do
diff --git a/electrons.irp.f b/electrons.irp.f
index c9fd51a..59e274f 100644
--- a/electrons.irp.f
+++ b/electrons.irp.f
@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ integer, nelec ]
  BEGIN_DOC
  ! Number of electrons
  END_DOC
- nelec = 2
+ nelec = 100
 END_PROVIDER
 
 
@@ -43,14 +43,14 @@ BEGIN_PROVIDER [double precision, factor_ee]
  BEGIN_DOC
  ! Electron-electron contribution to Jastrow factor
  END_DOC
- integer :: i, j
+ integer :: i, j, p
  double precision :: pow_ser = 0.0d0
  factor_ee = 0.0d0
 
  do j = 1 , nelec
     do i = 1, nelec
        do p = 2, nbord
-          pow_ser = pow_ser + bord_vect(p + 1) * rescale_ee(i, j) ** p
+          pow_ser = pow_ser + bord_vect(p) * rescale_ee(i, j) ** p
        end do
        factor_ee = factor_ee + bord_vect(1) * rescale_ee(i, j) &
             / (1 + bord_vect(2) * rescale_ee(i, j)) + pow_ser
diff --git a/jastrow b/jastrow
index dd3601f..353128b 100755
Binary files a/jastrow and b/jastrow differ
diff --git a/jastrow.irp.f b/jastrow.irp.f
index 3986c32..fd34fc7 100644
--- a/jastrow.irp.f
+++ b/jastrow.irp.f
@@ -1,6 +1,6 @@
 program jastrow
   implicit none
   print *, 'The total Jastrow factor'
-  print *, dexp(factor_ee + factor_en + factor_een)
+  print *, jastrow_full
 
 end program
diff --git a/jastrow_provider.irp.f b/jastrow_provider.irp.f
new file mode 100644
index 0000000..d127f57
--- /dev/null
+++ b/jastrow_provider.irp.f
@@ -0,0 +1,13 @@
+BEGIN_PROVIDER [ double precision, jastrow_full ]
+ implicit none
+ BEGIN_DOC
+ ! Complete jastrow factor 
+ END_DOC
+ print *, factor_ee
+ print *, factor_en
+ print *, factor_een
+
+ jastrow_full = dexp(factor_ee + factor_en + factor_een)
+
+END_PROVIDER
+
diff --git a/nuclei.irp.f b/nuclei.irp.f
index 568636b..c6352a3 100644
--- a/nuclei.irp.f
+++ b/nuclei.irp.f
@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ integer, nnuc ]
  BEGIN_DOC
  ! Number of nuclei
  END_DOC
- nnuc = 2
+ nnuc = 10
 END_PROVIDER
 
 
@@ -21,7 +21,7 @@ BEGIN_PROVIDER [ double precision, nuc_coord, (nnuc, 3) ]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, elnuc_dist, (nnuc, nnuc) ]
+BEGIN_PROVIDER [ double precision, elnuc_dist, (nelec, nnuc) ]
  implicit none
  BEGIN_DOC
  ! e-n distance
@@ -50,7 +50,7 @@ BEGIN_PROVIDER [double precision, factor_en]
  do j = 1 , nnuc
     do i = 1, nnuc
        do p = 2, naord
-          pow_ser = pow_ser + aord_vect(p + 1) * rescale_en(i, j) ** p
+          pow_ser = pow_ser + aord_vect(p) * rescale_en(i, j) ** p
        end do
        factor_en = factor_en + aord_vect(1) * rescale_en(i, j) &
             / (1 + aord_vect(2) * rescale_en(i, j)) + pow_ser
diff --git a/orders.irp.f b/orders.irp.f
index ab48763..d50dc44 100644
--- a/orders.irp.f
+++ b/orders.irp.f
@@ -2,24 +2,24 @@ BEGIN_PROVIDER [integer, naord]
  implicit none
  BEGIN_DOC
  ! Expansion order for f_en
- END_DOC 
- naord = 5   
+ END_DOC
+ naord = 5
 END_PROVIDER
 
 BEGIN_PROVIDER [integer, nbord]
  implicit none
  BEGIN_DOC
  ! Expansion order for f_ee
- END_DOC 
- nbord = 5   
+ END_DOC
+ nbord = 5
 END_PROVIDER
 
 BEGIN_PROVIDER [integer, ncord]
  implicit none
  BEGIN_DOC
  ! Expansion order for f_een
- END_DOC 
- ncord = 5   
+ END_DOC
+ ncord = 5
 END_PROVIDER
 
 BEGIN_PROVIDER [double precision, aord_vect, (naord)]
@@ -27,9 +27,9 @@ BEGIN_PROVIDER [double precision, aord_vect, (naord)]
  BEGIN_DOC
  ! Vector of the `a' coefficients
  END_DOC
- do i = 1, naord 
-    call random_number(aord_vect)
- end do
+ integer :: i
+ call random_number(aord_vect)
+ aord_vect = aord_vect*.1e-2
 END_PROVIDER
 
 BEGIN_PROVIDER [double precision, bord_vect, (nbord)]
@@ -37,17 +37,16 @@ BEGIN_PROVIDER [double precision, bord_vect, (nbord)]
  BEGIN_DOC
  ! Vector of the `b' coefficients
  END_DOC
- do i = 1, nbord 
-    call random_number(bord_vect)
- end do
+ integer :: i
+ call random_number(bord_vect)
+ bord_vect = bord_vect*.1e-6
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, cord_vect, (ncord)]
+BEGIN_PROVIDER [double precision, cord_vect, (0:ncord,0:ncord,ncord,nnuc)]
  implicit none
  BEGIN_DOC
  ! Vector of the `c' coefficients
  END_DOC
- do i = 1, ncord 
-    call random_number(cord_vect)
- end do
+ call random_number(cord_vect)
+ cord_vect = cord_vect*.1e-4
 END_PROVIDER
diff --git a/rescale.irp.f b/rescale.irp.f
index d919b59..69b7185 100644
--- a/rescale.irp.f
+++ b/rescale.irp.f
@@ -40,7 +40,7 @@ BEGIN_PROVIDER [ double precision, rescale_en, (nelec, nnuc) ]
  enddo
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, 3)]
+BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec)]
  implicit none
  BEGIN_DOC
  ! R = exp(-kappa r) for electron-electron for $J_{een}$
@@ -53,7 +53,7 @@ BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, 3)]
  enddo
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, rescale_een_n, (nnuc, 3)]
+BEGIN_PROVIDER [double precision, rescale_een_n, (nelec, nnuc)]
  implicit none
  BEGIN_DOC
  ! R = exp(-kappa r) for electron-electron for $J_{een}$