diff --git a/Makefile b/Makefile
index c94629d..9ae4c71 100644
--- a/Makefile
+++ b/Makefile
@@ -1,8 +1,8 @@
-IRPF90 = irpf90 -a -d
+IRPF90 = irpf90  #-a -d
 FC     = gfortran
 FCFLAGS= -O2 -ffree-line-length-none -I .
 NINJA  = ninja
-AR     = ar
+ARCHIVE= ar crs
 RANLIB = ranlib
 
 SRC=
diff --git a/el_nuc_el.irp.f b/el_nuc_el.irp.f
index 081693a..2f342ee 100644
--- a/el_nuc_el.irp.f
+++ b/el_nuc_el.irp.f
@@ -3,28 +3,60 @@ BEGIN_PROVIDER [double precision, factor_een]
  BEGIN_DOC
  ! Electron-electron nucleus contribution to Jastrow factor
  END_DOC
- integer :: i, j, alpha, p, k, l, lmax = 0
+ 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)
        do i = 1, nelec
-        do p = 2, ncord
-           do k = p - 1, 0
-              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)
-                 end if
-              end do
-           end do
-        end do
+          u = rescale_een_e(i, j)
+          a = rescale_een_n(i, alpha)
+          a2 = a * a
+          b2 = b * b
+          c = rescale_een_n(i, alpha) * rescale_een_n(j, alpha)
+          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
diff --git a/elec_coord.txt b/elec_coord.txt
new file mode 100644
index 0000000..50e632f
--- /dev/null
+++ b/elec_coord.txt
@@ -0,0 +1,10 @@
+ -0.250655104764153       0.503070975550133       -0.166554344502303     
+ -0.587812193472177      -0.128751981129274        0.187773606533075     
+  1.61335569047166       -0.615556732874863       -1.43165470979934     
+ -4.901239896295210E-003 -1.120440036458986E-002   1.99761909330422     
+  0.766647499681200      -0.293515395797937        3.66454589201239     
+ -0.127732483187947      -0.138975497694196       -8.669850480215846E-002
+ -0.232271834949124      -1.059321673434182E-002  -0.504862241464867     
+  1.09360863531826       -2.036103063808752E-003  -2.702796910818986E-002
+ -0.108090166832043       0.189161729653261        2.15398313919894     
+  0.397978144318712      -0.254277292595981        2.54553335476344     
diff --git a/electrons.irp.f b/electrons.irp.f
index c9fd51a..9ab69b5 100644
--- a/electrons.irp.f
+++ b/electrons.irp.f
@@ -3,7 +3,15 @@ BEGIN_PROVIDER [ integer, nelec ]
  BEGIN_DOC
  ! Number of electrons
  END_DOC
- nelec = 2
+ nelec = 10
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, nelec_up ]
+ implicit none
+ BEGIN_DOC
+ ! Number of alpha and beta electrons
+ END_DOC
+ nelec_up = 5
 END_PROVIDER
 
 
@@ -12,12 +20,16 @@ BEGIN_PROVIDER [ double precision, elec_coord, (nelec, 3) ]
  BEGIN_DOC
  ! Electron coordinates
  END_DOC
- integer :: i,j
- do j = 1 , 3
-   do i = 1, nelec
-     call random_number(elec_coord(i, j))
-   enddo
- enddo
+ character(len=*), parameter :: FILE_NAME = "elec_coord.txt"
+ integer :: fu, rc, i, j
+
+ open(action='read', file=FILE_NAME, iostat=rc, newunit=fu)
+
+ do i = 1, nelec
+    read(fu, *) elec_coord(i, :)
+ end do
+
+ close(fu)
 
 END_PROVIDER
 
@@ -26,15 +38,16 @@ BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ]
  BEGIN_DOC
  ! e-e distance
  END_DOC
- integer :: i,j
- double precision :: x,y,z
+ integer :: i, j
+ double precision :: x, y, z
+
  do j = 1, nelec
-  do i = 1, nelec
-    x = elec_coord(i, 1) - elec_coord(j, 1)
-    y = elec_coord(i, 2) - elec_coord(j, 2)
-    z = elec_coord(i, 3) - elec_coord(j, 3)
-    elec_dist(i, j) = dsqrt( x*x + y*y + z*z )
-  enddo
+    do i = 1, nelec
+       x = elec_coord(i, 1) - elec_coord(j, 1)
+       y = elec_coord(i, 2) - elec_coord(j, 2)
+       z = elec_coord(i, 3) - elec_coord(j, 3)
+       elec_dist(i, j) = dsqrt( x*x + y*y + z*z )
+    enddo
  enddo
 END_PROVIDER
 
@@ -43,18 +56,35 @@ BEGIN_PROVIDER [double precision, factor_ee]
  BEGIN_DOC
  ! Electron-electron contribution to Jastrow factor
  END_DOC
- integer :: i, j
- double precision :: pow_ser = 0.0d0
+ integer :: i, j, p, ipar
+ double precision :: pow_ser, x, spin_fact
+
  factor_ee = 0.0d0
 
- do j = 1 , nelec
+ do j = 1, nelec
     do i = 1, nelec
+       x = rescale_ee(i, j) 
+       pow_ser = 0.0d0
+       spin_fact = 1.0d0
+       ipar = 0
+
        do p = 2, nbord
-          pow_ser = pow_ser + bord_vect(p + 1) * rescale_ee(i, j) ** p
+          x = x * rescale_ee(i, j)
+          pow_ser = pow_ser + bord_vect(p + 1) * x
        end do
-       factor_ee = factor_ee + bord_vect(1) * rescale_ee(i, j) &
-            / (1 + bord_vect(2) * rescale_ee(i, j)) + pow_ser
+
+       if ((i.le.nelec_up .and. j.le.nelec_up) .or. &
+	       (i.gt.nelec_up .and. j.gt.nelec_up)) then
+           spin_fact = 0.5d0
+	   ipar = 1
+       end if
+
+       factor_ee = factor_ee + spin_fact * bord_vect(1) * rescale_ee(i, j) &
+            / (1.0d0 + bord_vect(2) * rescale_ee(i, j)) + pow_ser
+
     end do
  end do
+
  factor_ee = 0.5d0 * factor_ee
+
 END_PROVIDER
diff --git a/geometry.txt b/geometry.txt
new file mode 100644
index 0000000..9c8041c
--- /dev/null
+++ b/geometry.txt
@@ -0,0 +1,2 @@
+0.000000  0.000000  0.000000
+0.000000  0.000000  2.059801
diff --git a/jast_coeffs.txt b/jast_coeffs.txt
new file mode 100644
index 0000000..5b23c66
--- /dev/null
+++ b/jast_coeffs.txt
@@ -0,0 +1,35 @@
+0.00000000
+0.00000000
+-0.380512
+-0.157996
+-0.031558
+0.021512
+0.5000000
+0.153660
+0.0672262
+0.021570
+0.0073096
+0.002866
+0.571702
+-0.5142530
+-0.513043
+0.009486
+-0.004205
+0.4263258
+0.0828815
+0.0051186
+-0.0029978
+-0.0052704
+-0.000075
+-0.0830165
+0.0145434
+0.0514351
+0.000925
+-0.0040991
+0.0043276
+-0.00165447
+0.002614
+-0.001477
+-0.0011370
+-0.04010475
+0.00610671
diff --git a/jastrow b/jastrow
index dd3601f..c47738f 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..70a0acd
--- /dev/null
+++ b/jastrow_provider.irp.f
@@ -0,0 +1,15 @@
+BEGIN_PROVIDER [ double precision, jastrow_full ]
+ implicit none
+ BEGIN_DOC
+ ! Complete jastrow factor 
+ END_DOC
+ integer :: i, j
+
+ 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..3f3ebbb 100644
--- a/nuclei.irp.f
+++ b/nuclei.irp.f
@@ -7,21 +7,36 @@ BEGIN_PROVIDER [ integer, nnuc ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ integer, typenuc ]
+&BEGIN_PROVIDER [integer, typenuc_arr, (nnuc)]
+ implicit none
+ BEGIN_DOC
+ ! Number of nuclei
+ END_DOC
+ typenuc = 1
+ typenuc_arr = (/1, 1/)
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ double precision, nuc_coord, (nnuc, 3) ]
  implicit none
  BEGIN_DOC
  ! Nuclei coordinates
  END_DOC
- integer :: i, j
- do j = 1 , 3
-   do i = 1, nnuc
-     call random_number(nuc_coord(i, j))
-   enddo
- enddo
+ character(len=*), parameter :: FILE_NAME = "geometry.txt"
+ integer :: fu, rc, i
+ 
+ open(action='read', file=FILE_NAME, iostat=rc, newunit=fu)
+
+ do i = 1, nnuc
+    read(fu, *) nuc_coord(i, :)
+ end do
+ 
+ close(fu)
 
 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
@@ -29,12 +44,12 @@ BEGIN_PROVIDER [ double precision, elnuc_dist, (nnuc, nnuc) ]
  integer :: i, j
  double precision :: x, y, z
  do j = 1, nnuc
-  do i = 1, nelec
-    x = elec_coord(i, 1) - nuc_coord(j, 1)
-    y = elec_coord(i, 2) - nuc_coord(j, 2)
-    z = elec_coord(i, 3) - nuc_coord(j, 3)
-    elnuc_dist(i, j) = dsqrt( x*x + y*y + z*z )
-  enddo
+    do i = 1, nelec
+       x = elec_coord(i, 1) - nuc_coord(j, 1)
+       y = elec_coord(i, 2) - nuc_coord(j, 2)
+       z = elec_coord(i, 3) - nuc_coord(j, 3)
+       elnuc_dist(i, j) = dsqrt( x*x + y*y + z*z )
+    enddo
  enddo
 END_PROVIDER
 
@@ -43,19 +58,25 @@ BEGIN_PROVIDER [double precision, factor_en]
  BEGIN_DOC
  ! Electron-nuclei contribution to Jastrow factor
  END_DOC
- integer :: i, j, p
- double precision :: pow_ser = 0.0d0
+ integer :: i, j, p, q
+ double precision :: pow_ser, x
+
  factor_en = 0.0d0
 
  do j = 1 , nnuc
-    do i = 1, nnuc
+    do i = 1, nelec
+       x = rescale_en(i, j) 
+       pow_ser = 0.0d0
+
        do p = 2, naord
-          pow_ser = pow_ser + aord_vect(p + 1) * rescale_en(i, j) ** p
+          x = x * rescale_en(i, j) 
+          pow_ser = pow_ser + aord_vect(p + 1, typenuc_arr(j)) * x
        end do
-       factor_en = factor_en + aord_vect(1) * rescale_en(i, j) &
-            / (1 + aord_vect(2) * rescale_en(i, j)) + pow_ser
+
+       factor_en = factor_en + aord_vect(1, typenuc_arr(j)) * rescale_en(i, j) &
+            / (1 + aord_vect(2, typenuc_arr(j)) * rescale_en(i, j)) + pow_ser
+
     end do
  end do
 
- factor_en = 0.5d0 * factor_en
 END_PROVIDER
diff --git a/orders.irp.f b/orders.irp.f
index ab48763..187f1cb 100644
--- a/orders.irp.f
+++ b/orders.irp.f
@@ -2,52 +2,70 @@ 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)]
+ 
+BEGIN_PROVIDER [integer, dim_cord_vect]
  implicit none
  BEGIN_DOC
- ! Vector of the `a' coefficients
+ ! Recomputes the length of the unique C coefficients
  END_DOC
- do i = 1, naord 
-    call random_number(aord_vect)
- end do
-END_PROVIDER
+ integer :: k, p, l, lmax
 
-BEGIN_PROVIDER [double precision, bord_vect, (nbord)]
+ dim_cord_vect = 0
+
+ do p = 2, ncord
+    do k = 0, p - 1
+       if ( k /= 0 ) then
+          lmax = p - k
+       else
+          lmax = p - k - 2
+       end if
+       do l = iand(p - k, 1), lmax, 2
+          dim_cord_vect = dim_cord_vect + 1
+       end do
+    end do
+ end do
+ 
+END_PROVIDER
+ 
+
+BEGIN_PROVIDER [double precision, aord_vect, (naord + 1, typenuc)]
+&BEGIN_PROVIDER [double precision, bord_vect, (nbord + 1)]
+&BEGIN_PROVIDER [double precision, cord_vect, (dim_cord_vect, typenuc)]
  implicit none
  BEGIN_DOC
- ! Vector of the `b' coefficients
+ ! Read Jastow coefficients from file
  END_DOC
- do i = 1, nbord 
-    call random_number(bord_vect)
- end do
-END_PROVIDER
+ PROVIDE naord
+ PROVIDE nbord
+ PROVIDE ncord
+ character(len=*), parameter :: FILE_NAME = "jast_coeffs.txt"
+ integer :: i, fu, rc
+ 
+ open(action='read', file=FILE_NAME, iostat=rc, newunit=fu)
+ 
+ read(fu, *) aord_vect
+ read(fu, *) bord_vect
+ read(fu, *) cord_vect
+ 
+ close(fu)
 
-BEGIN_PROVIDER [double precision, cord_vect, (ncord)]
- implicit none
- BEGIN_DOC
- ! Vector of the `c' coefficients
- END_DOC
- do i = 1, ncord 
-    call random_number(cord_vect)
- end do
 END_PROVIDER
diff --git a/rescale.irp.f b/rescale.irp.f
index d919b59..56b0d66 100644
--- a/rescale.irp.f
+++ b/rescale.irp.f
@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, kappa ]
  BEGIN_DOC
  ! Constant in rescaling
  END_DOC
- kappa = 1.5d0
+ kappa = 0.6d0
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, kappa_inv ]
@@ -11,7 +11,7 @@ BEGIN_PROVIDER [ double precision, kappa_inv ]
  BEGIN_DOC
  ! inverse of kappa
  END_DOC
- kappa_inv = 1.d0 / kappa
+ kappa_inv = 1.0d0 / kappa
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, rescale_ee, (nelec, nelec) ]
@@ -20,9 +20,10 @@ BEGIN_PROVIDER [ double precision, rescale_ee, (nelec, nelec) ]
  ! R = (1 - exp(-kappa r))/kappa for electron-electron for $J_{ee}$
  END_DOC
  integer :: i, j
- do j=1,nelec
-   do i=1,nelec
-     rescale_ee(i, j) = (1.d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv
+
+ do j = 1, nelec
+   do i = 1, nelec
+     rescale_ee(i, j) = (1.0d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv
    enddo
  enddo
 END_PROVIDER
@@ -33,6 +34,7 @@ BEGIN_PROVIDER [ double precision, rescale_en, (nelec, nnuc) ]
  ! R = (1 - exp(-kappa r))/kappa for electron-nucleus for $J_{en}$
  END_DOC
  integer :: i, j
+
  do j = 1, nnuc
    do i = 1, nelec
      rescale_en(i, j) = (1.d0 - dexp(-kappa * elnuc_dist(i, j))) * kappa_inv
@@ -40,12 +42,13 @@ 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}$
  END_DOC
  integer :: i, j
+
  do j = 1, nelec
    do i = 1, nelec
      rescale_een_e(i, j) = dexp(-kappa * elec_dist(i, j))
@@ -53,12 +56,13 @@ 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}$
  END_DOC
  integer :: i, j
+
  do j = 1, nnuc
    do i = 1, nelec
       rescale_een_n(i, j) = dexp(-kappa * elnuc_dist(i, j))