Not working

src/BiInts/map_integrals.irp.f

@@ -32,21 +32,72 @@ end
 
 subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
   implicit none
-  integer, intent(out)           :: i,j,k,l
+  integer, intent(out)           :: i(8),j(8),k(8),l(8)
   integer*8, intent(in)          :: i1
   integer*8                      :: i2,i3
   real                           :: x
+  i = 0
   x = 0.5*(sqrt(8.*real(i1)+1.)-1.)
   i2 = ceiling(x)
   i3 = i1 - ishft(i2*i2-i2,-1)
 
-  l = 0.5*(sqrt(8.*real(i2)+1.)-1.)
-  l = ceiling(x)
-  j = i2 - ishft(l*l-l,-1)
+  x = 0.5*(sqrt(8.*real(i2)+1.)-1.)
+  l(1) = ceiling(x)
+  j(1) = i2 - ishft(l(1)*l(1)-l(1),-1)
 
   x = 0.5*(sqrt(8.*real(i3)+1.)-1.)
-  k = ceiling(x)
-  i = i3 - ishft(k*k-k,-1)
+  k(1) = ceiling(x)
+  i(1) = i3 - ishft(k(1)*k(1)-k(1),-1)
+
+              !ijkl
+  i(2) = k(1) !kjil
+  j(2) = j(1)
+  k(2) = i(1)
+  l(2) = l(1)
+
+  i(3) = i(1) !ilkj
+  j(3) = l(1)
+  k(3) = k(1)
+  l(3) = j(1)
+
+  i(4) = k(1) !klij
+  j(4) = l(1)
+  k(4) = i(1)
+  l(4) = j(1)
+
+  i(5) = j(1) !jilk
+  j(5) = i(1)
+  k(5) = l(1)
+  l(5) = k(1)
+
+  i(6) = l(1) !lijk
+  j(6) = i(1)
+  k(6) = j(1)
+  l(6) = k(1)
+
+  i(7) = j(1) !jkli
+  j(7) = k(1)
+  k(7) = l(1)
+  l(7) = i(1)
+
+  i(8) = l(1) !lkji
+  j(8) = k(1)
+  k(8) = j(1)
+  l(8) = i(1)
+
+  integer :: ii, jj
+  do ii=2,8
+    do jj=1,ii-1
+      if ( (i(ii) == i(jj)).and. &
+           (j(ii) == j(jj)).and. &
+           (k(ii) == k(jj)).and. &
+           (l(ii) == l(jj)) ) then
+         i(jj) = 0
+         exit
+      endif
+    enddo
+  enddo
+
 end
 
 

src/Hartree_Fock/Fock_matrix.irp.f

@@ -101,6 +101,7 @@ END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_alpha, (ao_num_align, ao_num) ]
 &BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_beta ,  (ao_num_align, ao_num) ]
+ use map_module
  implicit none
  BEGIN_DOC
  ! Alpha Fock matrix in AO basis set
@@ -159,7 +160,6 @@ END_PROVIDER
    !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,&
    !$OMP  ao_num_align,ao_bi_elec_integral_beta,HF_density_matrix_ao_alpha, &
    !$OMP  HF_density_matrix_ao_beta,lck) 
-   allocate(ao_ints_idx(ao_num_align),ao_ints_val(ao_num_align))
    !$OMP DO 
    do j=1,ao_num
      !DIR$ VECTOR ALIGNED
@@ -169,54 +169,42 @@ END_PROVIDER
      enddo
    enddo
    !$OMP END DO
-   !$OMP DO SCHEDULE(GUIDED)
-   do j=1,ao_num
-     do l=1,ao_num
-       do i=1,j-1
-         call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax)
-         call OMP_SET_LOCK(lck(i))
-         call OMP_SET_LOCK(lck(j))
-         do k1=1,kmax
-           k = ao_ints_idx(k1)
-           integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * ao_ints_val(k1)
-           ao_bi_elec_integral_alpha(i,j) += integral
-           ao_bi_elec_integral_beta (i,j) += integral
-           ao_bi_elec_integral_alpha(j,i) += integral
-           ao_bi_elec_integral_beta (j,i) += integral
-         enddo
-         do k1=1,kmax
-           k = ao_ints_idx(k1)
-           integral = ao_ints_val(k1)
-           ao_bi_elec_integral_alpha(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral
-           ao_bi_elec_integral_beta (l,j) -= HF_density_matrix_ao_beta (k,i) * integral
-           ao_bi_elec_integral_alpha(l,i) -= HF_density_matrix_ao_alpha(k,j) * integral
-           ao_bi_elec_integral_beta (l,i) -= HF_density_matrix_ao_beta (k,j) * integral
-         enddo
-         call OMP_UNSET_LOCK(lck(i))
-         call OMP_UNSET_LOCK(lck(j))
-       enddo
-       do i=j,j
-         call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax)
-         call OMP_SET_LOCK(lck(i))
-         do k1=1,kmax
-           k = ao_ints_idx(k1)
-           integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * ao_ints_val(k1)
-           ao_bi_elec_integral_alpha(i,j) += integral
-           ao_bi_elec_integral_beta (i,j) += integral
-         enddo
-         do k1=1,kmax
-           k = ao_ints_idx(k1)
-           integral = ao_ints_val(k1)
-           ao_bi_elec_integral_alpha(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral
-           ao_bi_elec_integral_beta (l,j) -= HF_density_matrix_ao_beta (k,i) * integral
-         enddo
-         call OMP_UNSET_LOCK(lck(i))
+   !$OMP END PARALLEL
+
+   integer(cache_map_size_kind)   :: n_elements_max, n_elements
+   integer(key_kind), allocatable :: keys(:)
+   double precision, allocatable  :: values(:)
+   call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max)
+   allocate(keys(n_elements_max), values(n_elements_max))
+
+   integer*8                      :: i8
+   integer                        :: ii(8), jj(8), kk(8), ll(8), k2
+
+   do i8=1,ao_integrals_map%map_size
+     n_elements = n_elements_max
+     call get_cache_map(ao_integrals_map,i8,keys,values,n_elements)
+     if (n_elements == 0) then
+       cycle
+     endif
+     do k1=1,n_elements
+       call bielec_integrals_index_reverse(kk,ii,ll,jj,keys(k1))
+       do k2=1,8
+         if (kk(k2) == 0) cycle
+         i = ii(k2)
+         j = jj(k2)
+         k = kk(k2)
+         l = ll(k2)
+         integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * values(k1)
+         ao_bi_elec_integral_alpha(i,j) += integral
+         ao_bi_elec_integral_beta (i,j) += integral
+!         integral = values(k1)
+!         ao_bi_elec_integral_alpha(i,l) -= HF_density_matrix_ao_alpha(k,l) * integral
+!         ao_bi_elec_integral_beta (i,l) -= HF_density_matrix_ao_beta (k,l) * integral
        enddo
      enddo
    enddo
-   !$OMP END DO
-   deallocate(ao_ints_val,ao_ints_idx)
-   !$OMP END PARALLEL
+   deallocate(keys,values)
+
    do i=1,ao_num
      call omp_destroy_lock(lck(i))
    enddo