working on mo 2e int framework

src/mo_two_e_ints/integrals_3_index.irp.f

@@ -25,3 +25,33 @@
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [complex*16, big_array_coulomb_integrals_periodic, (mo_num,mo_num, mo_num)]
+&BEGIN_PROVIDER [complex*16, big_array_exchange_integrals_periodic,(mo_num,mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! big_array_coulomb_integrals(j,i,k)  = <ij|kj> = (ik|jj)
+ ! big_array_exchange_integrals(j,i,k) = <ij|jk> = (ij|jk)
+ !    for both of these, i and k must be from same kpt for integral to be nonzero
+ ! TODO: only loop over half, and assign two elements:
+ !        b_a_coul_int(j,i,k) = b_a_coul_int(j,k,i)*
+ !        b_a_exch_int(j,i,k) = b_a_exch_int(j,k,i)*
+ END_DOC
+ integer :: i,j,k,l
+ complex*16 :: get_two_e_integral_periodic
+ complex*16 :: integral
+
+ do k = 1, mo_num
+  do i = 1, mo_num
+   do j = 1, mo_num
+     l = j
+     integral = get_two_e_integral_periodic(i,j,k,l,mo_integrals_map,mo_integrals_map_2)
+     big_array_coulomb_integrals(j,i,k) = integral
+     l = j
+     integral = get_two_e_integral_periodic(i,j,l,k,mo_integrals_map,mo_integrals_map_2)
+     big_array_exchange_integrals(j,i,k) = integral
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -28,9 +28,60 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
   integer(bit_kind)              :: mask_ijkl(N_int,4)
   integer(bit_kind)              :: mask_ijk(N_int,3)
   double precision               :: cpu_1, cpu_2, wall_1, wall_2
+  integer*8                      :: get_mo_map_size, mo_map_size
+  double precision, external     :: map_mb
 
   PROVIDE mo_class
 
+  if (is_periodic) then
+    mo_two_e_integrals_in_map = .True.
+    if (read_mo_two_e_integrals) then
+      print*,'Reading the MO integrals'
+      call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_periodic_1',mo_integrals_map)
+      call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_periodic_2',mo_integrals_map_2)
+      print*, 'MO integrals provided (periodic)'
+      return
+    else
+      PROVIDE ao_two_e_integrals_in_map
+    endif
+    
+    print *,  ''
+    print *,  'AO -> MO integrals transformation (periodic)'
+    print *,  '---------------------------------'
+    print *,  ''
+  
+    call wall_time(wall_1)
+    call cpu_time(cpu_1)
+  
+    if(no_vvvv_integrals)then
+      print*,'not implemented for periodic',irp_here
+      stop -1
+      call four_idx_novvvv_periodic
+    else
+      print*,'not implemented for periodic',irp_here
+      stop -1
+      call add_integrals_to_map_periodic(full_ijkl_bitmask_4)
+    endif
+  
+    call wall_time(wall_2)
+    call cpu_time(cpu_2)
+  
+    mo_map_size = get_mo_map_size()
+  
+    print*,'Molecular integrals provided:'
+    print*,' Size of MO map 1         ', map_mb(mo_integrals_map) ,'MB'
+    print*,' Size of MO map 2         ', map_mb(mo_integrals_map_2) ,'MB'
+    print*,' Number of MO integrals: ',  mo_map_size
+    print*,' cpu  time :',cpu_2 - cpu_1, 's'
+    print*,' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
+  
+    if (write_mo_two_e_integrals.and.mpi_master) then
+      call ezfio_set_work_empty(.False.)
+      call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_periodic_1',mo_integrals_map)
+      call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_periodic_2',mo_integrals_map_2)
+      call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
+    endif
+  else
   mo_two_e_integrals_in_map = .True.
   if (read_mo_two_e_integrals) then
     print*,'Reading the MO integrals'
@@ -58,10 +109,8 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
   call wall_time(wall_2)
   call cpu_time(cpu_2)
 
-  integer*8                      :: get_mo_map_size, mo_map_size
   mo_map_size = get_mo_map_size()
 
-  double precision, external     :: map_mb
   print*,'Molecular integrals provided:'
   print*,' Size of MO map           ', map_mb(mo_integrals_map) ,'MB'
   print*,' Number of MO integrals: ',  mo_map_size
@@ -73,6 +122,7 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
     call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
     call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
   endif
+  endif
 
 END_PROVIDER
 

src/utils_periodic/qp2-pbc-diff.txt

@@ -54,6 +54,18 @@ mo_two_e_ints
 	incomplete
 
 
+NOTES:
+  number of unique 4-tuples with 8-fold symmetry is a8(n)=n*(n+1)*(n^2+n+2)/8
+  number of unique 4-tuples with 4-fold symmetry is a4(n)=n^2*(n^2+3)/4
+  a8 is number of unique real 2e ints with n mos
+  a4 is number of unique* complex 2e ints with n mos (where p+i*q and p-i*q are counted as one, not two)
+  a4(n) = a8(n) + a8(n-1)
+  
+  we can already generate the list of <ij|kl> with unique values for the 8-fold case
+  the set of these for 4-fold symmetry is the union of the 8-fold set for n and the 8-fold set for n-1 with a simple transformation
+  <ij|kl>_{4,n} = <ij|kl>_{8,n} + <(k+1)j|i(l+1)>_{8,n-1}
+
+
 
 
 ############################