started working on complex mo 2e ints

src/mo_two_e_ints/map_integrals.irp.f

@@ -4,6 +4,7 @@ use map_module
 !! ======
 
 BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
+&BEGIN_PROVIDER [ type(map_type), mo_integrals_map_2 ]
   implicit none
   BEGIN_DOC
   ! MO integrals
@@ -11,9 +12,17 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
   integer(key_kind)              :: key_max
   integer(map_size_kind)         :: sze
   call two_e_integrals_index(mo_num,mo_num,mo_num,mo_num,key_max)
+  if (is_periodic) then
+    sze = key_max*2
+    call map_init(mo_integrals_map,sze)
+    call map_init(mo_integrals_map_2,sze)
+    print*,  'MO maps initialized (complex): ', 2*sze
+  else
     sze = key_max
     call map_init(mo_integrals_map,sze)
+    call map_init(mo_integrals_map_2,1_map_size_kind)
     print*, 'MO map initialized: ', sze
+  endif
 END_PROVIDER
 
 subroutine insert_into_mo_integrals_map(n_integrals,                 &
@@ -32,6 +41,22 @@ subroutine insert_into_mo_integrals_map(n_integrals,                 &
   call map_update(mo_integrals_map, buffer_i, buffer_values, n_integrals, thr)
 end
 
+subroutine insert_into_mo_integrals_map_2(n_integrals,               &
+      buffer_i, buffer_values, thr)
+  use map_module
+  implicit none
+
+  BEGIN_DOC
+  ! Create new entry into MO map, or accumulate in an existing entry
+  END_DOC
+
+  integer, intent(in)                :: n_integrals
+  integer(key_kind), intent(inout)   :: buffer_i(n_integrals)
+  real(integral_kind), intent(inout) :: buffer_values(n_integrals)
+  real(integral_kind), intent(in)    :: thr
+  call map_update(mo_integrals_map_2, buffer_i, buffer_values, n_integrals, thr)
+end
+
  BEGIN_PROVIDER [ integer*4, mo_integrals_cache_min ]
 &BEGIN_PROVIDER [ integer*4, mo_integrals_cache_max ]
 &BEGIN_PROVIDER [ integer*8, mo_integrals_cache_min_8 ]
@@ -85,6 +110,44 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*12
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ complex*16, mo_integrals_cache_periodic, (0_8:128_8*128_8*128_8*128_8) ]
+ implicit none
+ BEGIN_DOC
+ ! Cache of MO integrals for fast access
+ END_DOC
+ PROVIDE mo_two_e_integrals_in_map
+ integer*8                      :: i,j,k,l
+ integer*4                      :: i4,j4,k4,l4
+ integer*8                      :: ii
+ integer(key_kind)              :: idx
+ complex(integral_kind)            :: integral
+ complex*16 :: get_mo_two_e_integrals_periodic_simple
+ FREE ao_integrals_cache
+ !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral)
+ do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+   l4 = int(l,4)
+   do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+     k4 = int(k,4)
+     do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+       j4 = int(j,4)
+       do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+         i4 = int(i,4)
+         !DIR$ FORCEINLINE
+         integral = get_mo_two_e_integrals_periodic_simple(i,j,k,l,&
+                    mo_integrals_map,mo_integrals_map_2)
+         ii = l-mo_integrals_cache_min_8
+         ii = ior( shiftl(ii,7), k-mo_integrals_cache_min_8)
+         ii = ior( shiftl(ii,7), j-mo_integrals_cache_min_8)
+         ii = ior( shiftl(ii,7), i-mo_integrals_cache_min_8)
+         mo_integrals_cache_periodic(ii) = integral
+       enddo
+     enddo
+   enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+END_PROVIDER
+
 
 double precision function get_two_e_integral(i,j,k,l,map)
   use map_module
@@ -118,6 +181,39 @@ double precision function get_two_e_integral(i,j,k,l,map)
   endif
 end
 
+complex*16 function get_two_e_integral_periodic(i,j,k,l,map,map2)
+  use map_module
+  implicit none
+!  BEGIN_DOC
+!  ! Returns one integral <ij|kl> in the MO basis
+!  ! TODO: finish this
+!  END_DOC
+!  integer, intent(in)            :: i,j,k,l
+!  integer(key_kind)              :: idx
+!  integer                        :: ii
+!  integer*8                      :: ii_8
+!  type(map_type), intent(inout)  :: map
+!  real(integral_kind)            :: tmp
+!  PROVIDE mo_two_e_integrals_in_map mo_integrals_cache_periodic
+!  ii = l-mo_integrals_cache_min
+!  ii = ior(ii, k-mo_integrals_cache_min)
+!  ii = ior(ii, j-mo_integrals_cache_min)
+!  ii = ior(ii, i-mo_integrals_cache_min)
+!  if (iand(ii, -128) /= 0) then
+!    !DIR$ FORCEINLINE
+!    call two_e_integrals_index(i,j,k,l,idx)
+!    !DIR$ FORCEINLINE
+!    call map_get(map,idx,tmp)
+!    get_two_e_integral_periodic = dble(tmp)
+!  else
+!    ii_8 = int(l,8)-mo_integrals_cache_min_8
+!    ii_8 = ior( shiftl(ii_8,7), int(k,8)-mo_integrals_cache_min_8)
+!    ii_8 = ior( shiftl(ii_8,7), int(j,8)-mo_integrals_cache_min_8)
+!    ii_8 = ior( shiftl(ii_8,7), int(i,8)-mo_integrals_cache_min_8)
+!    get_two_e_integral = mo_integrals_cache(ii_8)
+!  endif
+end
+
 
 double precision function mo_two_e_integral(i,j,k,l)
   implicit none