Use flattened array

el_nuc_el.irp.f

@@ -3,55 +3,62 @@ 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
+ integer :: i, j, alpha, p, k, l, lmax, cidx
  double precision :: x, y, z, t, c_inv, u, a, b, a2, b2, c, t0
+
  PROVIDE cord_vect
  factor_een = 0.0d0
+ cidx = 0
 
  do alpha = 1, nnuc
     do j = 1, nelec
        b = rescale_een_n(j, alpha)
        do i = 1, nelec
-        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
-        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
-                ! factor_een = factor_een + cord_vect(l, k, p, alpha) * (y + z) * t
-                factor_een = factor_een + cord_vect(l + k + p + alpha) * (y + z) * t
-                t = t * c_inv
-                y = y * a2
-                z = z * b2
-              end do
-              x = x * u
-           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
+          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
+                   ! factor_een = factor_een + cord_vect(l, k, p, alpha) * (y + z) * t
+                   ! cidx = l + (ncord + 1) * k + (ncord + 1) * (ncord + 1) * p + &
+                   !      (ncord + 1) * (ncord + 1) * ncord * alpha
+                   ! here I try to use the flattened version of the array
+                   cidx = l + 6 * k + 6 * 6 * p + 6 * 6 * 5 * alpha
+                   print *, cidx
+                   factor_een = factor_een + cord_vect(cidx) * (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

jastrow_provider.irp.f

@@ -3,6 +3,18 @@ BEGIN_PROVIDER [ double precision, jastrow_full ]
  BEGIN_DOC
  ! Complete jastrow factor 
  END_DOC
+ ! integer :: i, j, k, l
+
+ ! do l = 1, nnuc
+ !    do k = 1, ncord
+ !       do j = 0, ncord
+ !          do i = 0, ncord
+ !             write(*, *) cord_vect_0(i, j, k, l)
+ !          end do
+ !       end do
+ !    end do
+ ! end do
+
 
  print *, factor_ee
  print *, factor_en

orders.irp.f

@@ -74,13 +74,13 @@ END_PROVIDER
 !  FREE seed
 ! END_PROVIDER
 ! 
-! BEGIN_PROVIDER [double precision, cord_vect, (0:ncord,0:ncord,ncord,nnuc)]
-!  implicit none
-!  BEGIN_DOC
-!  ! Vector of the `c' coefficients
-!  END_DOC
-!  PROVIDE seed
-!  call random_number(cord_vect)
-!  cord_vect = cord_vect*.1d-4
-!  FREE seed
-! END_PROVIDER
+BEGIN_PROVIDER [double precision, cord_vect_0, (0:ncord,0:ncord,ncord,nnuc)]
+ implicit none
+ BEGIN_DOC
+ ! Vector of the `c' coefficients
+ END_DOC
+ PROVIDE seed
+ call random_number(cord_vect_0)
+ cord_vect_0 = cord_vect_0 * .1d-4
+ FREE seed
+END_PROVIDER