From fee0ae8680fa36d376a69048e910acb9ebdd4c6c Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 2 Jun 2020 18:17:15 -0500 Subject: [PATCH] better mo transformation --- src/mo_two_e_ints/df_mo_ints.irp.f | 230 +++++++++++++++++++++++- src/mo_two_e_ints/mo_bi_integrals.irp.f | 3 +- 2 files changed, 228 insertions(+), 5 deletions(-) diff --git a/src/mo_two_e_ints/df_mo_ints.irp.f b/src/mo_two_e_ints/df_mo_ints.irp.f index 9bc1e3c0..3a61911e 100644 --- a/src/mo_two_e_ints/df_mo_ints.irp.f +++ b/src/mo_two_e_ints/df_mo_ints.irp.f @@ -19,6 +19,228 @@ 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 @@ -155,7 +377,7 @@ subroutine mo_map_fill_from_df_single 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_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_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 @@ -169,7 +391,7 @@ subroutine mo_map_fill_from_df_single 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_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_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 @@ -183,11 +405,11 @@ subroutine mo_map_fill_from_df_single if (n_integrals_1 > 0) then !call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1) - call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_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_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_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, & diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index a4b3530e..16322c19 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -44,7 +44,8 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] 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_single + !call mo_map_fill_from_df_single + call mo_map_fill_from_df_dot return else PROVIDE ao_two_e_integrals_in_map