Merge pull request #111 from kgasperich/features_kpts

Features kpts

src/mo_two_e_ints/df_mo_ints.irp.f

@@ -19,6 +19,441 @@ BEGIN_PROVIDER [complex*16, df_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_
 
 END_PROVIDER
 
+subroutine mo_map_fill_from_df_dot
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! fill mo bielec integral map using 3-index df integrals
+  END_DOC
+
+  integer :: i,k,j,l,mu
+  integer :: ki,kk,kj,kl
+  integer :: ii,ik,ij,il
+  integer :: kikk2,kjkl2,jl2,ik2
+  integer :: i_mo,j_mo,i_df
+
+  complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:)
+
+  complex*16 :: integral,mjl,mik
+  integer                        :: n_integrals_1, n_integrals_2
+  integer                        :: size_buffer
+  integer(key_kind),allocatable  :: buffer_i_1(:), buffer_i_2(:)
+  real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:)
+  double precision               :: tmp_re,tmp_im
+  integer                        :: mo_num_kpt_2
+
+  double precision               :: cpu_1, cpu_2, wall_1, wall_2, wall_0
+  double precision               :: map_mb
+
+  logical :: use_map1
+  integer(key_kind) :: idx_tmp
+  double precision :: sign
+  complex*16, external :: zdotc
+
+  mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt
+
+  size_buffer = min(mo_num_per_kpt*mo_num_per_kpt*mo_num_per_kpt,16000000)
+  print*, 'Providing the mo_bielec integrals from 3-index df integrals'
+  call write_time(6)
+!  call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write')
+!  TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals 
+
+  call wall_time(wall_1)
+  call cpu_time(cpu_1)
+
+  allocate( ints_jl(df_num,mo_num_per_kpt,mo_num_per_kpt))
+  allocate( ints_ik(df_num,mo_num_per_kpt,mo_num_per_kpt))
+
+  wall_0 = wall_1
+  do kl=1, kpt_num
+    do kj=1, kl
+      call idx2_tri_int(kj,kl,kjkl2)
+      if (kj < kl) then
+        do i_mo=1,mo_num_per_kpt
+          do j_mo=1,mo_num_per_kpt
+            do i_df=1,df_num
+              ints_jl(i_df,i_mo,j_mo) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kjkl2))
+            enddo
+          enddo
+        enddo
+      else
+        do i_mo=1,mo_num_per_kpt
+          do j_mo=1,mo_num_per_kpt
+            do i_df=1,df_num
+              ints_jl(i_df,i_mo,j_mo) = df_mo_integrals_complex(i_mo,j_mo,i_df,kjkl2)
+            enddo
+          enddo
+        enddo
+      endif
+
+      do kk=1,kl
+        ki=kconserv(kl,kk,kj)
+        if (ki>kl) cycle
+        call idx2_tri_int(ki,kk,kikk2)
+        if (ki < kk) then
+          do i_mo=1,mo_num_per_kpt
+            do j_mo=1,mo_num_per_kpt
+              do i_df=1,df_num
+                ints_ik(i_df,i_mo,j_mo) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kikk2))
+              enddo
+            enddo
+          enddo
+!          ints_ik = conjg(reshape(df_mo_integral_array(:,:,:,kikk2),(/mo_num_per_kpt,mo_num_per_kpt,df_num/),order=(/2,1,3/)))
+        else
+          do i_mo=1,mo_num_per_kpt
+            do j_mo=1,mo_num_per_kpt
+              do i_df=1,df_num
+                ints_ik(i_df,i_mo,j_mo) = df_mo_integrals_complex(i_mo,j_mo,i_df,kikk2)
+              enddo
+            enddo
+          enddo
+        endif
+
+        !$OMP PARALLEL PRIVATE(i,k,j,l,ii,ik,ij,il,jl2,ik2, &
+            !$OMP  mu, mik, mjl, &
+            !$OMP  n_integrals_1, buffer_i_1, buffer_values_1, &
+            !$OMP  n_integrals_2, buffer_i_2, buffer_values_2, &
+            !$OMP  idx_tmp, tmp_re, tmp_im, integral,sign,use_map1) &
+            !$OMP  DEFAULT(NONE)  &
+            !$OMP  SHARED(size_buffer, kpt_num, df_num, mo_num_per_kpt, mo_num_kpt_2, &
+            !$OMP  kl,kj,kjkl2,ints_jl, &
+            !$OMP  ki,kk,kikk2,ints_ik, &
+            !$OMP  kconserv, df_mo_integrals_complex, mo_integrals_threshold, &
+            !$OMP  mo_integrals_map, mo_integrals_map_2)
+
+        allocate( &
+          buffer_i_1(size_buffer), &
+          buffer_i_2(size_buffer), &
+          buffer_values_1(size_buffer), &
+          buffer_values_2(size_buffer) &
+        )
+      
+        n_integrals_1=0
+        n_integrals_2=0
+        !$OMP DO SCHEDULE(guided)
+        do il=1,mo_num_per_kpt
+          l=il+(kl-1)*mo_num_per_kpt
+          do ij=1,mo_num_per_kpt
+            j=ij+(kj-1)*mo_num_per_kpt
+            if (j>l) exit
+            call idx2_tri_int(j,l,jl2)
+            do ik=1,mo_num_per_kpt
+              k=ik+(kk-1)*mo_num_per_kpt
+              if (k>l) exit
+              do ii=1,mo_num_per_kpt
+                i=ii+(ki-1)*mo_num_per_kpt
+                if ((j==l) .and. (i>k)) exit
+                call idx2_tri_int(i,k,ik2)
+                if (ik2 > jl2) exit
+                integral = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
+!                print*,i,k,j,l,real(integral),imag(integral)
+                if (cdabs(integral) < mo_integrals_threshold) then
+                  cycle
+                endif
+                call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign)
+                tmp_re = dble(integral)
+                tmp_im = dimag(integral)
+                if (use_map1) then
+                  n_integrals_1 += 1
+                  buffer_i_1(n_integrals_1)=idx_tmp
+                  buffer_values_1(n_integrals_1)=tmp_re
+                  if (sign.ne.0.d0) then
+                    n_integrals_1 += 1
+                    buffer_i_1(n_integrals_1)=idx_tmp+1
+                    buffer_values_1(n_integrals_1)=tmp_im*sign
+                  endif 
+                  if (n_integrals_1 >= size(buffer_i_1)-1) then
+                    call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1)
+                    !call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold)
+                    n_integrals_1 = 0
+                  endif
+                else
+                  n_integrals_2 += 1
+                  buffer_i_2(n_integrals_2)=idx_tmp
+                  buffer_values_2(n_integrals_2)=tmp_re
+                  if (sign.ne.0.d0) then
+                    n_integrals_2 += 1
+                    buffer_i_2(n_integrals_2)=idx_tmp+1
+                    buffer_values_2(n_integrals_2)=tmp_im*sign
+                  endif
+                  if (n_integrals_2 >= size(buffer_i_2)-1) then
+                    call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2)
+                    !call insert_into_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold)
+                    n_integrals_2 = 0
+                  endif
+                endif
+
+              enddo !ii
+            enddo !ik
+          enddo !ij
+        enddo !il
+        !$OMP END DO NOWAIT
+
+        if (n_integrals_1 > 0) then
+          call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1)
+          !call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold)
+        endif
+        if (n_integrals_2 > 0) then
+          call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2)
+          !call insert_into_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold)
+        endif
+        deallocate( &
+          buffer_i_1, &
+          buffer_i_2, &
+          buffer_values_1, &
+          buffer_values_2 &
+          )
+        !$OMP END PARALLEL
+      enddo !kk
+    enddo !kj
+    call wall_time(wall_2)
+    if (wall_2 - wall_0 > 1.d0) then
+      wall_0 = wall_2
+      print*, 100.*float(kl)/float(kpt_num), '% in ',               &
+          wall_2-wall_1,'s',map_mb(mo_integrals_map),'+',map_mb(mo_integrals_map_2),'MB'
+    endif
+
+  enddo !kl
+  deallocate( ints_jl,ints_ik )
+
+  call map_sort(mo_integrals_map)
+  call map_unique(mo_integrals_map)
+  call map_sort(mo_integrals_map_2)
+  call map_unique(mo_integrals_map_2)
+  !call map_merge(mo_integrals_map)
+  !call map_merge(mo_integrals_map_2)
+
+  !!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_1',mo_integrals_map)
+  !!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_2',mo_integrals_map_2)
+  !!call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
+  
+  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()
+  
+  print*,'MO integrals provided:'
+  print*,' Size of MO map           ', map_mb(mo_integrals_map),'+',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), ')'
+  
+end subroutine mo_map_fill_from_df_dot
+
+subroutine mo_map_fill_from_df_single
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! fill mo bielec integral map using 3-index df integrals
+  END_DOC
+
+  integer :: i,k,j,l,mu
+  integer :: ki,kk,kj,kl
+  integer :: ii,ik,ij,il
+  integer :: kikk2,kjkl2,jl2,ik2
+  integer :: i_mo,j_mo,i_df
+
+  complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:)
+
+  complex*16 :: integral,mjl,mik
+  integer                        :: n_integrals_1, n_integrals_2
+  integer                        :: size_buffer
+  integer(key_kind),allocatable  :: buffer_i_1(:), buffer_i_2(:)
+  real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:)
+  double precision               :: tmp_re,tmp_im
+  integer                        :: mo_num_kpt_2
+
+  double precision               :: cpu_1, cpu_2, wall_1, wall_2, wall_0
+  double precision               :: map_mb
+
+  logical :: use_map1
+  integer(key_kind) :: idx_tmp
+  double precision :: sign
+
+  mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt
+
+  size_buffer = min(mo_num_per_kpt*mo_num_per_kpt*mo_num_per_kpt,16000000)
+  print*, 'Providing the mo_bielec integrals from 3-index df integrals'
+  call write_time(6)
+!  call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write')
+!  TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals 
+
+  call wall_time(wall_1)
+  call cpu_time(cpu_1)
+
+  allocate( ints_jl(mo_num_per_kpt,mo_num_per_kpt,df_num))
+  allocate( ints_ik(mo_num_per_kpt,mo_num_per_kpt,df_num))
+
+  wall_0 = wall_1
+  do kl=1, kpt_num
+    do kj=1, kl
+      call idx2_tri_int(kj,kl,kjkl2)
+      if (kj < kl) then
+        do i_mo=1,mo_num_per_kpt
+          do j_mo=1,mo_num_per_kpt
+            do i_df=1,df_num
+              ints_jl(i_mo,j_mo,i_df) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kjkl2))
+            enddo
+          enddo
+        enddo
+      else
+        ints_jl = df_mo_integrals_complex(:,:,:,kjkl2)
+      endif
+
+      do kk=1,kl
+        ki=kconserv(kl,kk,kj)
+        if (ki>kl) cycle
+        call idx2_tri_int(ki,kk,kikk2)
+        if (ki < kk) then
+          do i_mo=1,mo_num_per_kpt
+            do j_mo=1,mo_num_per_kpt
+              do i_df=1,df_num
+                ints_ik(i_mo,j_mo,i_df) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kikk2))
+              enddo
+            enddo
+          enddo
+!          ints_ik = conjg(reshape(df_mo_integral_array(:,:,:,kikk2),(/mo_num_per_kpt,mo_num_per_kpt,df_num/),order=(/2,1,3/)))
+        else
+          ints_ik = df_mo_integrals_complex(:,:,:,kikk2)
+        endif
+
+        !$OMP PARALLEL PRIVATE(i,k,j,l,ii,ik,ij,il,jl2,ik2, &
+            !$OMP  mu, mik, mjl, &
+            !$OMP  n_integrals_1, buffer_i_1, buffer_values_1, &
+            !$OMP  n_integrals_2, buffer_i_2, buffer_values_2, &
+            !$OMP  idx_tmp, tmp_re, tmp_im, integral,sign,use_map1) &
+            !$OMP  DEFAULT(NONE)  &
+            !$OMP  SHARED(size_buffer, kpt_num, df_num, mo_num_per_kpt, mo_num_kpt_2, &
+            !$OMP  kl,kj,kjkl2,ints_jl, &
+            !$OMP  ki,kk,kikk2,ints_ik, &
+            !$OMP  kconserv, df_mo_integrals_complex, mo_integrals_threshold, &
+            !$OMP  mo_integrals_map, mo_integrals_map_2)
+
+        allocate( &
+          buffer_i_1(size_buffer), &
+          buffer_i_2(size_buffer), &
+          buffer_values_1(size_buffer), &
+          buffer_values_2(size_buffer) &
+        )
+      
+        n_integrals_1=0
+        n_integrals_2=0
+        !$OMP DO SCHEDULE(guided)
+        do mu=1,df_num
+          do il=1,mo_num_per_kpt
+            l=il+(kl-1)*mo_num_per_kpt
+            do ij=1,mo_num_per_kpt
+              j=ij+(kj-1)*mo_num_per_kpt
+              if (j>l) exit
+              call idx2_tri_int(j,l,jl2)
+              mjl = ints_jl(ij,il,mu)
+              if (mjl.eq.(0.d0,0.d0)) cycle
+              do ik=1,mo_num_per_kpt
+                k=ik+(kk-1)*mo_num_per_kpt
+                if (k>l) exit
+                do ii=1,mo_num_per_kpt
+                  i=ii+(ki-1)*mo_num_per_kpt
+                  if ((j==l) .and. (i>k)) exit
+                  call idx2_tri_int(i,k,ik2)
+                  if (ik2 > jl2) exit
+                  mik = ints_ik(ii,ik,mu)
+                  integral = mik * dconjg(mjl)
+!                  print*,i,k,j,l,real(integral),imag(integral)
+                  if (cdabs(integral) < mo_integrals_threshold) then
+                    cycle
+                  endif
+                  call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign)
+                  tmp_re = dble(integral)
+                  tmp_im = dimag(integral)
+                  if (use_map1) then
+                    n_integrals_1 += 1
+                    buffer_i_1(n_integrals_1)=idx_tmp
+                    buffer_values_1(n_integrals_1)=tmp_re
+                    if (sign.ne.0.d0) then
+                      n_integrals_1 += 1
+                      buffer_i_1(n_integrals_1)=idx_tmp+1
+                      buffer_values_1(n_integrals_1)=tmp_im*sign
+                    endif 
+                    if (n_integrals_1 >= size(buffer_i_1)-1) then
+                      !call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1)
+                      call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold)
+                      n_integrals_1 = 0
+                    endif
+                  else
+                    n_integrals_2 += 1
+                    buffer_i_2(n_integrals_2)=idx_tmp
+                    buffer_values_2(n_integrals_2)=tmp_re
+                    if (sign.ne.0.d0) then
+                      n_integrals_2 += 1
+                      buffer_i_2(n_integrals_2)=idx_tmp+1
+                      buffer_values_2(n_integrals_2)=tmp_im*sign
+                    endif
+                    if (n_integrals_2 >= size(buffer_i_2)-1) then
+                      !call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2)
+                      call insert_into_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold)
+                      n_integrals_2 = 0
+                    endif
+                  endif
+
+                enddo !ii
+              enddo !ik
+            enddo !ij
+          enddo !il
+        enddo !mu
+        !$OMP END DO NOWAIT
+
+        if (n_integrals_1 > 0) then
+          !call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1)
+          call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold)
+        endif
+        if (n_integrals_2 > 0) then
+          !call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2)
+          call insert_into_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold)
+        endif
+        deallocate( &
+          buffer_i_1, &
+          buffer_i_2, &
+          buffer_values_1, &
+          buffer_values_2 &
+          )
+        !$OMP END PARALLEL
+      enddo !kk
+    enddo !kj
+    call wall_time(wall_2)
+    if (wall_2 - wall_0 > 1.d0) then
+      wall_0 = wall_2
+      print*, 100.*float(kl)/float(kpt_num), '% in ',               &
+          wall_2-wall_1,'s',map_mb(mo_integrals_map),'+',map_mb(mo_integrals_map_2),'MB'
+    endif
+
+  enddo !kl
+  deallocate( ints_jl,ints_ik )
+
+  !call map_sort(mo_integrals_map)
+  !call map_unique(mo_integrals_map)
+  !call map_sort(mo_integrals_map_2)
+  !call map_unique(mo_integrals_map_2)
+  call map_merge(mo_integrals_map)
+  call map_merge(mo_integrals_map_2)
+  !!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_1',mo_integrals_map)
+  !!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_2',mo_integrals_map_2)
+  !!call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
+  
+  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()
+  
+  print*,'MO integrals provided:'
+  print*,' Size of MO map           ', map_mb(mo_integrals_map),'+',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), ')'
+  
+end subroutine mo_map_fill_from_df_single
+
 subroutine mo_map_fill_from_df
   use map_module
   implicit none

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -43,7 +43,9 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
       return
     else if (read_df_mo_integrals.or.read_df_ao_integrals) then
       PROVIDE df_mo_integrals_complex
-      call mo_map_fill_from_df
+      !call mo_map_fill_from_df
+      !call mo_map_fill_from_df_single
+      call mo_map_fill_from_df_dot
       return
     else
       PROVIDE ao_two_e_integrals_in_map

src/utils_complex/export_integrals_ao_cplx.irp.f

@@ -159,7 +159,7 @@ provide ao_two_e_integrals_in_map
    if (cdabs(tmp_cmplx-int2e_tmp1).gt.thr0) then
      print'(4(I4),4(E15.7))',i,j,k,l,tmp_cmplx,int2e_tmp1
    endif
-    integer*8 :: ii
+    integer :: ii
           ii = l-ao_integrals_cache_min
           ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
           ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)