Optimizations: reduce number of exponentials

Makefile

@@ -1,8 +1,9 @@
-IRPF90 = irpf90  #-a -d
+IRPF90 = irpf90 --codelet=factor_een:100000 #-a -d
-FC     = gfortran
+FC     = ifort -xHost -g
 FCFLAGS= -O2 -ffree-line-length-none -I .
 NINJA  = ninja
 AR = ar 
+ARCHIVE = ar crs
 RANLIB = ranlib
 
 SRC=

el_nuc_el.irp.f

@@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, factor_een ]
  ! ElectronE-electron-nuclei contribution to Jastrow factor
  END_DOC
  integer :: i, j, a, p, k, l, lmax, m
- double precision :: riam, rjam_cn, rial, rjal, rijk
+ double precision :: rjam_cn
  double precision :: cn
 
  factor_een = 0.0d0
@@ -21,14 +21,16 @@ BEGIN_PROVIDER [ double precision, factor_een ]
           m = (p - k - l) / 2
           do a = 1, nnuc
              cn = cord_vect_lkp(l, k, p, typenuc_arr(a))
-             do j = 1, nelec
+             rjam_cn = rescale_een_n(2, a, m) * cn
-                rjal = rescale_een_n(j, a, l)
+             factor_een = factor_een + rescale_een_e(1,2,k) * &
+              (rescale_een_n(1,a,l) + rescale_een_n(2,a,l)) * &
+               rescale_een_n(1,a,m) * rjam_cn
+             do j = 3, nelec
                 rjam_cn = rescale_een_n(j, a, m) * cn
                 do i = 1, j - 1
-                   rial = rescale_een_n(i, a, l)
+                   factor_een = factor_een + rescale_een_e(i,j,k) * &
-                   riam = rescale_een_n(i, a, m)
+                    (rescale_een_n(i,a,l) + rescale_een_n(j,a,l)) * &
-                   rijk = rescale_een_e(i, j, k)
+                     rescale_een_n(i,a,m) * rjam_cn
-                   factor_een = factor_een + rijk * (rial + rjal) * riam * rjam_cn
                 enddo
              enddo
           enddo

rescale.irp.f

@@ -20,11 +20,15 @@ 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
+ double precision :: x
 
  do j = 1, nelec
-    do i = 1, nelec
+    do i = 1, j-1
-       rescale_ee(i, j) = (1.0d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv
+       x = (1.0d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv
+       rescale_ee(i, j) = x
+       rescale_ee(j, i) = x
     enddo
+    rescale_ee(j, j) = 0.d0
  enddo
 END_PROVIDER
 
@@ -36,9 +40,11 @@ BEGIN_PROVIDER [ double precision, rescale_ee_deriv_e, (4, nelec, nelec) ]
  ! Dimension 4 : d2x + d2y + d2z
  END_DOC
  integer :: i, j, ii
+ double precision :: f
 
  do j = 1, nelec
     do i = 1, nelec
+       f = 1.d0 - kappa*rescale_ee(i,j) ! == dexp(-kappa * elec_dist(i, j))
        do ii = 1, 4
           rescale_ee_deriv_e(ii, i, j) = elec_dist_deriv_e(ii, i, j)
        end do
@@ -48,7 +54,7 @@ BEGIN_PROVIDER [ double precision, rescale_ee_deriv_e, (4, nelec, nelec) ]
        (-kappa * rescale_ee_deriv_e(3, i, j) * rescale_ee_deriv_e(3, i, j))
        do ii = 1, 4
           rescale_ee_deriv_e(ii, i, j) = rescale_ee_deriv_e(ii, i, j) &
-               * dexp(-kappa * elec_dist(i, j))
+               * f
        enddo
     enddo
  enddo
@@ -76,9 +82,11 @@ BEGIN_PROVIDER [ double precision, rescale_en_deriv_e, (4, nelec, nnuc) ]
  ! Dimension 4 : d2x + d2y + d2z
  END_DOC
  integer :: i, ii, a
+ double precision :: f
 
  do a = 1, nnuc
     do i = 1, nelec
+       f = 1.d0 - kappa*rescale_en(i,a) ! == dexp(-kappa * elnuc_dist(i, a))
        do ii = 1, 4
           rescale_en_deriv_e(ii, i, a) = elnuc_dist_deriv_e(ii, i, a)
        end do
@@ -88,7 +96,7 @@ BEGIN_PROVIDER [ double precision, rescale_en_deriv_e, (4, nelec, nnuc) ]
        (-kappa * rescale_en_deriv_e(3, i, a) * rescale_en_deriv_e(3, i, a))
        do ii = 1, 4
           rescale_en_deriv_e(ii, i, a) = rescale_en_deriv_e(ii, i, a) &
-               * dexp(-kappa * elnuc_dist(i, a))
+               * f
        enddo
     enddo
  enddo
@@ -100,18 +108,28 @@ BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec, 0:ncord)]
  ! R = exp(-kappa r) for electron-electron for $J_{een}$
  END_DOC
  integer :: i, j, l
- double precision :: kappa_l
+ double precision :: x
+ double precision, parameter :: f = dexp(1.d0)
+
+ rescale_een_e(:, :, 0) = 1.d0
 
- do l = 0, ncord
-    kappa_l = -dble(l) * kappa
  do j = 1, nelec
-       do i = 1, nelec
+    do i = 1, j-1
-          rescale_een_e(i, j, l) = kappa_l * elec_dist(i, j)
+       x = dexp(-kappa * elec_dist(i, j))
-       enddo
+       rescale_een_e(i, j, 1) = x
+       rescale_een_e(j, i, 1) = x
     enddo
  enddo
 
- rescale_een_e = dexp(rescale_een_e)
+ do l = 2, ncord
+    do j = 1, nelec
+       do i = 1, j-1
+          x = rescale_een_e(i, j, l-1) * rescale_een_e(i, j, 1)
+          rescale_een_e(i, j, l) = x
+          rescale_een_e(j, i, l) = x
+       enddo
+    enddo
+ enddo
 
  do l = 0, ncord
     do j = 1, nelec
@@ -126,18 +144,24 @@ BEGIN_PROVIDER [double precision, rescale_een_n, (nelec, nnuc, 0:ncord)]
  ! R = exp(-kappa r) for electron-electron for $J_{een}$
  END_DOC
  integer :: i, a, l
- double precision :: kappa_l
+ double precision :: x
+ double precision, parameter :: f = dexp(1.d0)
+
+ rescale_een_n(:,:,0) = 1.d0
 
- do l = 0, ncord
-    kappa_l = - dble(l) * kappa
  do a = 1, nnuc
     do i = 1, nelec
-          rescale_een_n(i, a, l) = kappa_l * elnuc_dist(i, a)
+       rescale_een_n(i, a, 1) = dexp(-kappa * elnuc_dist(i, a))
-       enddo
     enddo
  enddo
 
- rescale_een_n = dexp(rescale_een_n)
+ do l = 2, ncord
+    do a = 1, nnuc
+       do i = 1, nelec
+          rescale_een_n(i, a, l) = rescale_een_n(i, a, l-1) * rescale_een_n(i, a, 1)
+       enddo
+    enddo
+ enddo
 
 END_PROVIDER
 
@@ -207,6 +231,7 @@ BEGIN_PROVIDER [double precision, rescale_een_e_deriv_e, (4, nelec, nelec, 0:nco
  integer :: i, ii, j, l
  double precision :: kappa_l
 
+ !TODO: Check if rescale_een_e_deriv_e(:,:,0) = 0.d0
  do l = 0, ncord
     kappa_l = - dble(l) * kappa
     do j = 1, nelec