From 4e933906323ef96abf570f594403b45fcf9ec611 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 22 Jan 2020 11:35:41 -0600 Subject: [PATCH] working on two e ints --- src/ao_one_e_ints/ao_overlap.irp.f | 11 +- src/ao_one_e_ints/kin_ao_ints.irp.f | 4 +- src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f | 30 ++ src/ao_two_e_ints/map_integrals.irp.f | 479 ++++++++---------- src/ao_two_e_ints/two_e_integrals.irp.f | 142 +++--- src/utils/map_module.f90 | 19 +- .../export_integrals_ao_periodic.irp.f | 212 ++++++++ .../import_integrals_ao_periodic.irp.f | 97 ++-- 8 files changed, 618 insertions(+), 376 deletions(-) create mode 100644 src/utils_periodic/export_integrals_ao_periodic.irp.f diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index d7300936..6510fd23 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -75,7 +75,16 @@ BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ] BEGIN_DOC ! Imaginary part of the overlap END_DOC - ao_overlap_imag = 0.d0 + if (read_ao_integrals_overlap) then + call ezfio_get_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num)) + print *, 'AO overlap integrals read from disk' + else + ao_overlap_imag = 0.d0 + endif + if (write_ao_integrals_overlap) then + call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num)) + print *, 'AO overlap integrals written to disk' + endif END_PROVIDER BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ] diff --git a/src/ao_one_e_ints/kin_ao_ints.irp.f b/src/ao_one_e_ints/kin_ao_ints.irp.f index 4f117deb..442c1f88 100644 --- a/src/ao_one_e_ints/kin_ao_ints.irp.f +++ b/src/ao_one_e_ints/kin_ao_ints.irp.f @@ -160,13 +160,13 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)] integer :: i,j,k,l if (read_ao_integrals_kinetic) then - call ezfio_get_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag) + call ezfio_get_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag) print *, 'AO kinetic integrals read from disk' else print *, irp_here, ': Not yet implemented' endif if (write_ao_integrals_kinetic) then - call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag) + call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag) print *, 'AO kinetic integrals written to disk' endif END_PROVIDER diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f index 988bbe0a..a92ba1f4 100644 --- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -27,6 +27,36 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)] END_PROVIDER +BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_imag, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Imaginary part of the pseudo_integrals + END_DOC + if (read_ao_integrals_pseudo) then + call ezfio_get_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num)) + print *, 'AO pseudo_integrals integrals read from disk' + else + ao_pseudo_integrals_imag = 0.d0 + endif + if (write_ao_integrals_pseudo) then + call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num)) + print *, 'AO pseudo_integrals integrals written to disk' + endif +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, ao_pseudo_integrals_complex, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Overlap for complex AOs + END_DOC + integer :: i,j + do j=1,ao_num + do i=1,ao_num + ao_pseudo_integrals_complex(i,j) = dcmplx( ao_pseudo_integrals(i,j), ao_pseudo_integrals_imag(i,j) ) + enddo + enddo +END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] implicit none BEGIN_DOC diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index e1a49357..e9d2740c 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -4,6 +4,7 @@ use map_module !! ====== BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] +&BEGIN_PROVIDER [ type(map_type), ao_integrals_map_2 ] implicit none BEGIN_DOC ! AO integrals @@ -11,22 +12,17 @@ BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] integer(key_kind) :: key_max integer(map_size_kind) :: sze call two_e_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max) - sze = key_max - call map_init(ao_integrals_map,sze) - print*, 'AO map initialized : ', sze -END_PROVIDER - -BEGIN_PROVIDER [ type(map_type), ao_integrals_map_periodic ] - implicit none - BEGIN_DOC - ! AO integrals - END_DOC - integer(key_kind) :: key_max - integer(map_size_kind) :: sze - call two_e_integrals_index_2fold(ao_num,ao_num,ao_num,ao_num,key_max) - sze = key_max - call map_init(ao_integrals_map_periodic,sze) - print*, 'complex AO map initialized : ', sze + if (is_periodic) then + sze = key_max*2 + call map_init(ao_integrals_map,sze) + call map_init(ao_integrals_map_2,sze) + print*, 'AO maps initialized (complex): ', 2*sze + else + sze = key_max + call map_init(ao_integrals_map,sze) + call map_init(ao_integrals_map_2,1_map_size_kind) + print*, 'AO map initialized : ', sze + endif END_PROVIDER subroutine two_e_integrals_index(i,j,k,l,i1) @@ -34,7 +30,7 @@ subroutine two_e_integrals_index(i,j,k,l,i1) implicit none BEGIN_DOC ! Gives a unique index for i,j,k,l using permtuation symmetry. -! i <-> k, j <-> l, and (i,k) <-> (j,l) for non-periodic systems +! i <-> k, j <-> l, and (i,k) <-> (j,l) END_DOC integer, intent(in) :: i,j,k,l integer(key_kind), intent(out) :: i1 @@ -50,6 +46,28 @@ subroutine two_e_integrals_index(i,j,k,l,i1) i1 = i1+shiftr(i2*i2-i2,1) end +subroutine two_e_integrals_index_periodic(i,j,k,l,i1,p,q) + use map_module + implicit none + BEGIN_DOC +! Gives a unique index for i,j,k,l using permtuation symmetry. +! i <-> k, j <-> l, and (i,k) <-> (j,l) + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer(key_kind) :: r,s,i2 + integer(key_kind),intent(out) :: p,q + p = min(i,k) + r = max(i,k) + p = p+shiftr(r*r-r,1) + q = min(j,l) + s = max(j,l) + q = q+shiftr(s*s-s,1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+shiftr(i2*i2-i2,1) +end + subroutine two_e_integrals_index_reverse(i,j,k,l,i1) @@ -143,152 +161,6 @@ end -subroutine ao_idx2_sq(i,j,ij) - implicit none - integer, intent(in) :: i,j - integer, intent(out) :: ij - if (ij) then - ij=(i-1)*(i-1)+2*j-mod(i,2) - else - ij=i*i - endif -end - -subroutine idx2_tri_int(i,j,ij) - implicit none - integer, intent(in) :: i,j - integer, intent(out) :: ij - integer :: p,q - p = max(i,j) - q = min(i,j) - ij = q+ishft(p*p-p,-1) -end - -subroutine ao_idx2_tri_key(i,j,ij) - use map_module - implicit none - integer, intent(in) :: i,j - integer(key_kind), intent(out) :: ij - integer(key_kind) :: p,q - p = max(i,j) - q = min(i,j) - ij = q+ishft(p*p-p,-1) -end - -subroutine two_e_integrals_index_2fold(i,j,k,l,i1) - use map_module - implicit none - integer, intent(in) :: i,j,k,l - integer(key_kind), intent(out) :: i1 - integer :: ik,jl - - call ao_idx2_sq(i,k,ik) - call ao_idx2_sq(j,l,jl) - call ao_idx2_tri_key(ik,jl,i1) -end - -subroutine ao_idx2_sq_rev(i,k,ik) - BEGIN_DOC - ! reverse square compound index - END_DOC -! p = ceiling(dsqrt(dble(ik))) -! q = ceiling(0.5d0*(dble(ik)-dble((p-1)*(p-1)))) -! if (mod(ik,2)==0) then -! k=p -! i=q -! else -! i=p -! k=q -! endif - integer, intent(in) :: ik - integer, intent(out) :: i,k - integer :: pq(0:1),i1,i2 - pq(0) = ceiling(dsqrt(dble(ik))) - pq(1) = ceiling(0.5d0*(dble(ik)-dble((pq(0)-1)*(pq(0)-1)))) - i1=mod(ik,2) - i2=mod(ik+1,2) - - k=pq(i1) - i=pq(i2) -end - -subroutine ao_idx2_tri_rev_key(i,k,ik) - use map_module - BEGIN_DOC - !return i<=k - END_DOC - integer(key_kind), intent(in) :: ik - integer, intent(out) :: i,k - integer(key_kind) :: tmp_k - k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) - tmp_k = k - i = int(ik - ishft(tmp_k*tmp_k-tmp_k,-1)) -end - -subroutine idx2_tri_rev_int(i,k,ik) - BEGIN_DOC - !return i<=k - END_DOC - integer, intent(in) :: ik - integer, intent(out) :: i,k - k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) - i = int(ik - ishft(k*k-k,-1)) -end - -subroutine two_e_integrals_index_reverse_2fold(i,j,k,l,i1) - use map_module - implicit none - integer, intent(out) :: i(2),j(2),k(2),l(2) - integer(key_kind), intent(in) :: i1 - integer(key_kind) :: i0 - integer :: i2,i3 - i = 0 - call ao_idx2_tri_rev_key(i3,i2,i1) - - call ao_idx2_sq_rev(j(1),l(1),i2) - call ao_idx2_sq_rev(i(1),k(1),i3) - - !ijkl - i(2) = j(1) !jilk - j(2) = i(1) - k(2) = l(1) - l(2) = k(1) - -! i(3) = k(1) !klij complex conjugate -! j(3) = l(1) -! k(3) = i(1) -! l(3) = j(1) -! -! i(4) = l(1) !lkji complex conjugate -! j(4) = k(1) -! k(4) = j(1) -! l(4) = i(1) - - integer :: ii - if ( (i(1)==i(2)).and. & - (j(1)==j(2)).and. & - (k(1)==k(2)).and. & - (l(1)==l(2)) ) then - i(2) = 0 - endif -! This has been tested with up to 1000 AOs, and all the reverse indices are -! correct ! We can remove the test -! do ii=1,2 -! if (i(ii) /= 0) then -! call two_e_integrals_index_2fold(i(ii),j(ii),k(ii),l(ii),i0) -! if (i1 /= i0) then -! print *, i1, i0 -! print *, i(ii), j(ii), k(ii), l(ii) -! stop 'two_e_integrals_index_reverse_2fold failed' -! endif -! endif -! enddo -end - - - BEGIN_PROVIDER [ integer, ao_integrals_cache_min ] &BEGIN_PROVIDER [ integer, ao_integrals_cache_max ] @@ -384,32 +256,47 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] real(integral_kind) :: tmp_re, tmp_im integer(key_kind) :: idx_re,idx_im complex(integral_kind) :: integral + integer(key_kind) :: p,q,r,s,ik,jl + logical :: ilek, jlel, iklejl - !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx1,idx2,tmp_re,tmp_im,idx_re,idx_im,ii,integral) + !$OMP PARALLEL DO PRIVATE (ilek,jlel,p,q,r,s, ik,jl,iklejl, & + !$OMP i,j,k,l,idx1,idx2,tmp_re,tmp_im,idx_re,idx_im,ii,integral) do l=ao_integrals_cache_min,ao_integrals_cache_max do k=ao_integrals_cache_min,ao_integrals_cache_max do j=ao_integrals_cache_min,ao_integrals_cache_max do i=ao_integrals_cache_min,ao_integrals_cache_max !DIR$ FORCEINLINE - call two_e_integrals_index_2fold(i,j,k,l,idx1) - !DIR$ FORCEINLINE - call two_e_integrals_index_2fold(k,l,i,j,idx2) - idx_re = min(idx1,idx2) - idx_im = max(idx1,idx2) - !DIR$ FORCEINLINE - call map_get(ao_integrals_map,idx_re,tmp_re) - if (idx_re /= idx_im) then - call map_get(ao_integrals_map,idx_im,tmp_im) - if (idx1 < idx2) then + call two_e_integrals_index(i,j,k,l,idx1) + ilek = (i.le.k) + jlel = (j.le.l) + idx1 = 2*idx1 - 1 + if (ilek.eqv.jlel) then !map1 + !TODO: merge these calls using map_get_2 + call map_get(ao_integrals_map,idx1,tmp_re) + call map_get(ao_integrals_map,idx1+1,tmp_im) + if (ilek) then integral = dcmplx(tmp_re,tmp_im) else integral = dcmplx(tmp_re,-tmp_im) endif - else - tmp_im = 0.d0 - integral = dcmplx(tmp_re,tmp_im) - endif + else !map2 + !TODO: merge these calls using map_get_2 + call map_get(ao_integrals_map_2,idx1,tmp_re) + call map_get(ao_integrals_map_2,idx1+1,tmp_im) + p = min(i,k) + r = max(i,k) + ik = p+shiftr(r*r-r,1) + q = min(j,l) + s = max(j,l) + jl = q+shiftr(s*s-s,1) + iklejl = (ik.le.jl) + if (ilek.eqv.iklejl) then + integral = dcmplx(tmp_re,tmp_im) + else + integral = dcmplx(tmp_re,-tmp_im) + endif + endif ii = l-ao_integrals_cache_min ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) @@ -424,8 +311,53 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] END_PROVIDER +subroutine ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx,sign) + use map_module + implicit none + BEGIN_DOC + ! get position of periodic AO integral + ! use_map1: true if integral is in first ao map, false if integral is in second ao map + ! idx: position of real part of integral in map (imag part is at idx+1) + ! sign: sign of imaginary part + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: idx + logical, intent(out) :: use_map1 + double precision, intent(out) :: sign + integer(key_kind) :: p,q,r,s,ik,jl + logical :: ilek, jlel, iklejl, ikeqjl + ! i.le.k, j.le.l, tri(i,k).le.tri(j,l) + !DIR$ FORCEINLINE + call two_e_integrals_index_periodic(i,j,k,l,idx,ik,jl) + ilek = (i.le.k) + jlel = (j.le.l) + idx = 2*idx-1 + ikeqjl = (ik.eq.jl) + if (ilek.eqv.jlel) then !map1 + use_map1=.True. + if (ikeqjl) then + sign=0.d0 + else if (ilek) then + sign=1.d0 + else + sign=-1.d0 + endif + else !map2 + use_map1=.False. + if (ikeqjl) then + sign=0.d0 + else + iklejl = (ik.le.jl) + if (ilek.eqv.iklejl) then + sign=1.d0 + else + sign=-1.d0 + endif + endif + endif +end -complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result) +complex*16 function get_ao_two_e_integral_periodic_simple(i,j,k,l,map,map2) result(result) use map_module implicit none BEGIN_DOC @@ -435,10 +367,66 @@ complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result) integer(key_kind) :: idx1,idx2 real(integral_kind) :: tmp_re, tmp_im integer(key_kind) :: idx_re,idx_im - type(map_type), intent(inout) :: map + type(map_type), intent(inout) :: map,map2 integer :: ii complex(integral_kind) :: tmp - PROVIDE ao_two_e_integrals_in_map_periodic ao_integrals_cache_periodic ao_integrals_cache_min + integer(key_kind) :: p,q,r,s,ik,jl + logical :: ilek, jlel, iklejl + ! a.le.c, b.le.d, tri(a,c).le.tri(b,d) + PROVIDE ao_two_e_integrals_in_map + !DIR$ FORCEINLINE + call two_e_integrals_index(i,j,k,l,idx1) + ilek = (i.le.k) + jlel = (j.le.l) + idx1 = idx1*2-1 + if (ilek.eqv.jlel) then !map1 + !TODO: merge these calls using map_get_2 + call map_get(map,idx1,tmp_re) + call map_get(map,idx1+1,tmp_im) + if (ilek) then + tmp = dcmplx(tmp_re,tmp_im) + else + tmp = dcmplx(tmp_re,-tmp_im) + endif + else !map2 + !TODO: merge these calls using map_get_2 + call map_get(map2,idx1,tmp_re) + call map_get(map2,idx1+1,tmp_im) + p = min(i,k) + r = max(i,k) + ik = p+shiftr(r*r-r,1) + q = min(j,l) + s = max(j,l) + jl = q+shiftr(s*s-s,1) + iklejl = (ik.le.jl) + if (ilek.eqv.iklejl) then + tmp = dcmplx(tmp_re,tmp_im) + else + tmp = dcmplx(tmp_re,-tmp_im) + endif + endif + result = tmp +end + + +complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map,map2) result(result) + use map_module + implicit none + BEGIN_DOC + ! Gets one AO bi-electronic integral from the AO map + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx1,idx2 + real(integral_kind) :: tmp_re, tmp_im + integer(key_kind) :: idx_re,idx_im + type(map_type), intent(inout) :: map,map2 + integer :: ii + complex(integral_kind) :: tmp + complex(integral_kind) :: get_ao_two_e_integral_periodic_simple + integer(key_kind) :: p,q,r,s,ik,jl + logical :: ilek, jlel, iklejl + ! a.le.c, b.le.d, tri(a,c).le.tri(b,d) + PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min !DIR$ FORCEINLINE if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then tmp = (0.d0,0.d0) @@ -450,25 +438,7 @@ complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result) ii = ior(ii, j-ao_integrals_cache_min) ii = ior(ii, i-ao_integrals_cache_min) if (iand(ii, -64) /= 0) then - !DIR$ FORCEINLINE - call two_e_integrals_index_2fold(i,j,k,l,idx1) - !DIR$ FORCEINLINE - call two_e_integrals_index_2fold(k,l,i,j,idx2) - idx_re = min(idx1,idx2) - idx_im = max(idx1,idx2) - !DIR$ FORCEINLINE - call map_get(ao_integrals_map,idx_re,tmp_re) - if (idx_re /= idx_im) then - call map_get(ao_integrals_map,idx_im,tmp_im) - if (idx1 < idx2) then - tmp = dcmplx(tmp_re,tmp_im) - else - tmp = dcmplx(tmp_re,-tmp_im) - endif - else - tmp_im = 0.d0 - tmp = dcmplx(tmp_re,tmp_im) - endif + tmp = get_ao_two_e_integral_periodic_simple(i,j,k,l,map,map2) else ii = l-ao_integrals_cache_min ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) @@ -528,14 +498,14 @@ subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val) PROVIDE ao_two_e_integrals_in_map ao_integrals_map thresh = ao_integrals_threshold - if (ao_overlap_abs_periodic(j,l) < thresh) then + if (ao_overlap_abs(j,l) < thresh) then out_val = (0.d0,0.d0) return endif complex*16 :: get_ao_two_e_integral_periodic do i=1,sze - out_val(i) = get_ao_two_e_integral_periodic(i,j,k,l,ao_integrals_map) + out_val(i) = get_ao_two_e_integral_periodic(i,j,k,l,ao_integrals_map,ao_integrals_map_2) enddo end @@ -554,6 +524,10 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z integer :: i integer(key_kind) :: hash double precision :: thresh,tmp + if(is_periodic) then + print*,'not implemented for periodic:',irp_here + stop -1 + endif PROVIDE ao_two_e_integrals_in_map thresh = ao_integrals_threshold @@ -598,6 +572,10 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out integer(key_kind) :: hash double precision :: tmp + if(is_periodic) then + print*,'not implemented for periodic:',irp_here + stop -1 + endif PROVIDE ao_two_e_integrals_in_map non_zero_int = 0 if (ao_overlap_abs(j,l) < thresh) then @@ -644,6 +622,10 @@ subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,s integer(key_kind) :: hash double precision :: tmp + if(is_periodic) then + print*,'not implemented for periodic:',irp_here + stop -1 + endif PROVIDE ao_two_e_integrals_in_map non_zero_int = 0 if (ao_overlap_abs(j,l) < thresh) then @@ -682,7 +664,7 @@ function get_ao_map_size() BEGIN_DOC ! Returns the number of elements in the AO map END_DOC - get_ao_map_size = ao_integrals_map % n_elements + get_ao_map_size = ao_integrals_map % n_elements + ao_integrals_map_2 % n_elements end subroutine clear_ao_map @@ -692,6 +674,9 @@ subroutine clear_ao_map END_DOC call map_deinit(ao_integrals_map) FREE ao_integrals_map + call map_deinit(ao_integrals_map_2) + FREE ao_integrals_map_2 + end @@ -709,82 +694,18 @@ subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values) call map_append(ao_integrals_map, buffer_i, buffer_values, n_integrals) end +subroutine insert_into_ao_integrals_map_2(n_integrals,buffer_i, buffer_values) + use map_module + implicit none + BEGIN_DOC + ! Create new entry into AO map + END_DOC -!subroutine dump_ao_integrals(filename) -! use map_module -! implicit none -! BEGIN_DOC -! ! Save to disk the |AO| integrals -! END_DOC -! character*(*), intent(in) :: filename -! integer(cache_key_kind), pointer :: key(:) -! real(integral_kind), pointer :: val(:) -! integer*8 :: i,j, n -! if (.not.mpi_master) then -! return -! endif -! call ezfio_set_work_empty(.False.) -! open(unit=66,file=filename,FORM='unformatted') -! write(66) integral_kind, key_kind -! write(66) ao_integrals_map%sorted, ao_integrals_map%map_size, & -! ao_integrals_map%n_elements -! do i=0_8,ao_integrals_map%map_size -! write(66) ao_integrals_map%map(i)%sorted, ao_integrals_map%map(i)%map_size,& -! ao_integrals_map%map(i)%n_elements -! enddo -! do i=0_8,ao_integrals_map%map_size -! key => ao_integrals_map%map(i)%key -! val => ao_integrals_map%map(i)%value -! n = ao_integrals_map%map(i)%n_elements -! write(66) (key(j), j=1,n), (val(j), j=1,n) -! enddo -! close(66) -! -!end + integer, intent(in) :: n_integrals + integer(key_kind), intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + + call map_append(ao_integrals_map_2, buffer_i, buffer_values, n_integrals) +end -!integer function load_ao_integrals(filename) -! implicit none -! BEGIN_DOC -! ! Read from disk the |AO| integrals -! END_DOC -! character*(*), intent(in) :: filename -! integer*8 :: i -! integer(cache_key_kind), pointer :: key(:) -! real(integral_kind), pointer :: val(:) -! integer :: iknd, kknd -! integer*8 :: n, j -! load_ao_integrals = 1 -! open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') -! read(66,err=98,end=98) iknd, kknd -! if (iknd /= integral_kind) then -! print *, 'Wrong integrals kind in file :', iknd -! stop 1 -! endif -! if (kknd /= key_kind) then -! print *, 'Wrong key kind in file :', kknd -! stop 1 -! endif -! read(66,err=98,end=98) ao_integrals_map%sorted, ao_integrals_map%map_size,& -! ao_integrals_map%n_elements -! do i=0_8, ao_integrals_map%map_size -! read(66,err=99,end=99) ao_integrals_map%map(i)%sorted, & -! ao_integrals_map%map(i)%map_size, ao_integrals_map%map(i)%n_elements -! call cache_map_reallocate(ao_integrals_map%map(i),ao_integrals_map%map(i)%map_size) -! enddo -! do i=0_8, ao_integrals_map%map_size -! key => ao_integrals_map%map(i)%key -! val => ao_integrals_map%map(i)%value -! n = ao_integrals_map%map(i)%n_elements -! read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) -! enddo -! call map_sort(ao_integrals_map) -! load_ao_integrals = 0 -! return -! 99 continue -! call map_deinit(ao_integrals_map) -! 98 continue -! stop 'Problem reading ao_integrals_map file in work/' -! -!end -! diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index a2bde897..e3ca0566 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -348,77 +348,91 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] integer :: kk, m, j1, i1, lmax character*(64) :: fmt - integral = ao_two_e_integral(1,1,1,1) double precision :: map_mb PROVIDE read_ao_two_e_integrals io_ao_two_e_integrals - if (read_ao_two_e_integrals) then - print*,'Reading the AO integrals' - call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) - print*, 'AO integrals provided' - ao_two_e_integrals_in_map = .True. - return - endif - - print*, 'Providing the AO integrals' - call wall_time(wall_0) - call wall_time(wall_1) - call cpu_time(cpu_1) - - integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull - call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals') - - character(len=:), allocatable :: task - allocate(character(len=ao_num*12) :: task) - write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' - do l=1,ao_num - write(task,fmt) (i,l, i=1,l) - integer, external :: add_task_to_taskserver - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then - stop 'Unable to add task to server' + if (is_periodic) then + if (read_ao_two_e_integrals) then + print*,'Reading the AO integrals (periodic)' + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_1',ao_integrals_map) + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_2',ao_integrals_map_2) + print*, 'AO integrals provided (periodic)' + ao_two_e_integrals_in_map = .True. + return + else + print*,'calculation of periodic AOs not implemented' + stop -1 endif - enddo - deallocate(task) - integer, external :: zmq_set_running - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif - - PROVIDE nproc - !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1) - i = omp_get_thread_num() - if (i==0) then - call ao_two_e_integrals_in_map_collector(zmq_socket_pull) - else - call ao_two_e_integrals_in_map_slave_inproc(i) + else + if (read_ao_two_e_integrals) then + print*,'Reading the AO integrals' + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + print*, 'AO integrals provided' + ao_two_e_integrals_in_map = .True. + return + endif + + integral = ao_two_e_integral(1,1,1,1) + print*, 'Providing the AO integrals' + call wall_time(wall_0) + call wall_time(wall_1) + call cpu_time(cpu_1) + + integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull + call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals') + + character(len=:), allocatable :: task + allocate(character(len=ao_num*12) :: task) + write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' + do l=1,ao_num + write(task,fmt) (i,l, i=1,l) + integer, external :: add_task_to_taskserver + if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then + stop 'Unable to add task to server' endif - !$OMP END PARALLEL - - call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals') - - - print*, 'Sorting the map' - call map_sort(ao_integrals_map) - call cpu_time(cpu_2) - call wall_time(wall_2) - integer(map_size_kind) :: get_ao_map_size, ao_map_size - ao_map_size = get_ao_map_size() - - print*, 'AO integrals provided:' - print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' - print*, ' Number of AO integrals :', ao_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+tiny(1.d0)), ' )' - - ao_two_e_integrals_in_map = .True. - - if (write_ao_two_e_integrals.and.mpi_master) then - call ezfio_set_work_empty(.False.) - call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) - call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + enddo + deallocate(task) + + integer, external :: zmq_set_running + if (zmq_set_running(zmq_to_qp_run_socket) == -1) then + print *, irp_here, ': Failed in zmq_set_running' + endif + + PROVIDE nproc + !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call ao_two_e_integrals_in_map_collector(zmq_socket_pull) + else + call ao_two_e_integrals_in_map_slave_inproc(i) + endif + !$OMP END PARALLEL + + call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals') + + + print*, 'Sorting the map' + call map_sort(ao_integrals_map) + call cpu_time(cpu_2) + call wall_time(wall_2) + integer(map_size_kind) :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + + print*, 'AO integrals provided:' + print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' + print*, ' Number of AO integrals :', ao_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+tiny(1.d0)), ' )' + + ao_two_e_integrals_in_map = .True. + + if (write_ao_two_e_integrals.and.mpi_master) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + endif endif - END_PROVIDER BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ] diff --git a/src/utils/map_module.f90 b/src/utils/map_module.f90 index 98e73470..e41ee8d0 100644 --- a/src/utils/map_module.f90 +++ b/src/utils/map_module.f90 @@ -531,13 +531,30 @@ subroutine map_get(map, key, value) real(integral_kind), intent(out) :: value integer(map_size_kind) :: idx_cache integer(cache_map_size_kind) :: idx - + idx=1 ! index in tha pointers array idx_cache = shiftr(key,map_shift) !DIR$ FORCEINLINE call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx) end +subroutine map_get_2(map, key, value1, value2) + use map_module + implicit none + type (map_type), intent(inout) :: map + integer(key_kind), intent(in) :: key + real(integral_kind), intent(out) :: value1, value2 + integer(map_size_kind) :: idx_cache + integer(cache_map_size_kind) :: idx + + idx=1 + ! index in tha pointers array + idx_cache = shiftr(key,map_shift) + !DIR$ FORCEINLINE + call cache_map_get_interval(map%map(idx_cache), key, value1, 1, map%map(idx_cache)%n_elements,idx) + call cache_map_get_interval(map%map(idx_cache), key+1, value2, idx+1, idx+2, idx) +end + subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx) use map_module implicit none diff --git a/src/utils_periodic/export_integrals_ao_periodic.irp.f b/src/utils_periodic/export_integrals_ao_periodic.irp.f new file mode 100644 index 00000000..8d60ff49 --- /dev/null +++ b/src/utils_periodic/export_integrals_ao_periodic.irp.f @@ -0,0 +1,212 @@ +program print_integrals + call run +end + +subroutine run + use map_module + implicit none + + integer :: iunit + integer :: getunitandopen + + integer ::i,j,k,l + double precision :: integral + double precision, allocatable :: A(:,:), B(:,:) + double precision :: tmp_re, tmp_im + + integer :: n_integrals_1, n_integrals_2 + integer(key_kind), allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind), allocatable :: buffer_values_1(:), buffer_values_2(:) + logical :: use_map1 + integer(key_kind) :: idx_tmp + double precision :: sign + + + +provide ao_two_e_integrals_in_map + allocate (A(ao_num,ao_num), B(ao_num,ao_num) ) + + A(1,1) = huge(1.d0) + iunit = getunitandopen('E.qp','r') + read (iunit,*,end=9) A(1,1) + 9 continue + close(iunit) + if (A(1,1) /= huge(1.d0)) then +! call ezfio_set_nuclei_nuclear_repulsion(A(1,1)) +! call ezfio_set_nuclei_io_nuclear_repulsion("Read") + print*, nuclear_repulsion,A(1,1) + endif + + A = 0.d0 + B = 0.d0 + ! iunit = getunitandopen('T.qp','r') + ! do + ! read (iunit,*,end=10) i,j, tmp_re, tmp_im + ! A(i,j) = tmp_re + ! B(i,j) = tmp_im + ! print*,ao_kinetic_integrals(i,j),A(i,j) + ! print*,ao_kinetic_integrals_imag(i,j),B(i,j) + ! if (i.ne.j) then + ! A(j,i) = tmp_re + ! B(j,i) = -tmp_im + ! print*,ao_kinetic_integrals(j,i),A(j,i) + ! print*,ao_kinetic_integrals_imag(j,i),B(j,i) + ! endif + ! enddo + ! 10 continue + ! close(iunit) +! call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A(1:ao_num, 1:ao_num)) +! call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(B(1:ao_num, 1:ao_num)) +! call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic("Read") + + A = 0.d0 + B = 0.d0 + ! iunit = getunitandopen('S.qp','r') + ! do + ! read (iunit,*,end=11) i,j, tmp_re, tmp_im + ! A(i,j) = tmp_re + ! B(i,j) = tmp_im + ! print*,real(ao_overlap_complex(i,j)),A(i,j) + ! print*,imag(ao_overlap_complex(i,j)),B(i,j) + ! print*,ao_overlap_imag(i,j),B(i,j) + ! if (i.ne.j) then + ! A(j,i) = tmp_re + ! B(j,i) = -tmp_im + ! print*,real(ao_overlap_complex(j,i)),A(j,i) + ! print*,imag(ao_overlap_complex(j,i)),B(j,i) + ! print*,ao_overlap_imag(j,i),B(j,i) + ! endif + ! enddo + ! 11 continue + ! close(iunit) +! call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A(1:ao_num, 1:ao_num)) +! call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(B(1:ao_num, 1:ao_num)) +! call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap("Read") + + A = 0.d0 + B = 0.d0 +! iunit = getunitandopen('P.qp','r') +! do +! read (iunit,*,end=14) i,j, tmp_re, tmp_im +! A(i,j) = tmp_re +! B(i,j) = tmp_im +! print*,ao_pseudo_integrals(i,j),A(i,j) +! print*,ao_pseudo_integrals_imag(i,j),B(i,j) +! ! print*,real(ao_integrals_pseudo(i,j)),A(i,j) +! ! print*,imag(ao_integrals_pseudo(i,j)),B(i,j) +! if (i.ne.j) then +! A(j,i) = tmp_re +! B(j,i) = -tmp_im +! print*,ao_pseudo_integrals(j,i),A(j,i) +! print*,ao_pseudo_integrals_imag(j,i),B(j,i) +! ! print*,real(ao_integrals_pseudo(j,i)),A(j,i) +! ! print*,imag(ao_integrals_pseudo(j,i)),B(j,i) +! endif +! enddo +! 14 continue +! close(iunit) +! call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A(1:ao_num,1:ao_num)) +! call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(B(1:ao_num,1:ao_num)) +! call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo("Read") + + A = 0.d0 + B = 0.d0 +! iunit = getunitandopen('V.qp','r') +! do +! read (iunit,*,end=12) i,j, tmp_re, tmp_im +! A(i,j) = tmp_re +! B(i,j) = tmp_im +! print*,ao_integrals_n_e(i,j),A(i,j) +! print*,ao_integrals_n_e_imag(i,j),B(i,j) +! if (i.ne.j) then +! A(j,i) = tmp_re +! B(j,i) = -tmp_im +! print*,ao_integrals_n_e(j,i),A(j,i) +! print*,ao_integrals_n_e_imag(j,i),B(j,i) +! endif +! enddo +! 12 continue +! close(iunit) +! call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A(1:ao_num, 1:ao_num)) +! call ezfio_set_ao_one_e_ints_ao_integrals_n_e_imag(B(1:ao_num, 1:ao_num)) +! call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e("Read") + complex*16 :: int2e_tmp1,int2e_tmp2,get_ao_two_e_integral_periodic_simple,get_ao_two_e_integral_periodic + double precision :: tmp3,tmp4,tmp5,tmp6 + allocate(buffer_i_1(ao_num**3), buffer_values_1(ao_num**3)) + allocate(buffer_i_2(ao_num**3), buffer_values_2(ao_num**3)) + iunit = getunitandopen('W.qp','r') + n_integrals_1=0 + n_integrals_2=0 + buffer_values_1 = 0.d0 + buffer_values_2 = 0.d0 + do + read (iunit,*,end=13) i,j,k,l, tmp_re, tmp_im + int2e_tmp1 = get_ao_two_e_integral_periodic_simple(i,j,k,l,ao_integrals_map,ao_integrals_map_2) + int2e_tmp2 = get_ao_two_e_integral_periodic(i,j,k,l,ao_integrals_map,ao_integrals_map_2) + print'(4(I4),3(E15.7))',i,j,k,l,tmp_re,real(int2e_tmp1),real(int2e_tmp2) + print'(4(I4),3(E15.7))',i,j,k,l,tmp_im,imag(int2e_tmp1),imag(int2e_tmp2) + call ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) + print*,use_map1,idx_tmp,sign + call map_get(ao_integrals_map,idx_tmp,tmp3) + call map_get(ao_integrals_map_2,idx_tmp,tmp4) + call map_get(ao_integrals_map,idx_tmp+1,tmp5) + call map_get(ao_integrals_map_2,idx_tmp+1,tmp6) + print*,tmp3,tmp4 + print*,tmp5,tmp6 +! if (use_map1) then +! n_integrals_1 += 1 +! buffer_i_1(n_integrals_1-1)=idx_tmp +! buffer_values_1(n_integrals_1-1)=tmp_re +! if (sign.ne.0.d0) then +! n_integrals_1 += 1 +! buffer_i_1(n_integrals_2)=idx_tmp+1 +! buffer_values_1(n_integrals_1)=tmp_im*sign +! endif +! if (n_integrals_1 >= size(buffer_i_1)-1) then +!! call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) +! n_integrals_1 = 0 +! endif +! else +! n_integrals_2 += 1 +! buffer_i_2(n_integrals_2-1)=idx_tmp +! buffer_values_2(n_integrals_2-1)=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 insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) +! n_integrals_2 = 0 +! endif +! endif + enddo + 13 continue + close(iunit) + +! if (n_integrals_1 > 0) then +!! call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) +! endif +! if (n_integrals_2 > 0) then +!! call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) +! endif +! +! call map_sort(ao_integrals_map) +! call map_unique(ao_integrals_map) +! call map_sort(ao_integrals_map_2) +! call map_unique(ao_integrals_map_2) +! +! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_1',ao_integrals_map) +! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_2',ao_integrals_map_2) +! call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read' +print*,'map1' + do i=0,ao_integrals_map%map_size + print*,i,ao_integrals_map%map(i)%value(:) + print*,i,ao_integrals_map%map(i)%key(:) + enddo +print*,'map2' + do i=0,ao_integrals_map_2%map_size + print*,i,ao_integrals_map_2%map(i)%value(:) + print*,i,ao_integrals_map_2%map(i)%key(:) + enddo +end diff --git a/src/utils_periodic/import_integrals_ao_periodic.irp.f b/src/utils_periodic/import_integrals_ao_periodic.irp.f index 79eb8fe0..039af101 100644 --- a/src/utils_periodic/import_integrals_ao_periodic.irp.f +++ b/src/utils_periodic/import_integrals_ao_periodic.irp.f @@ -17,9 +17,13 @@ subroutine run double precision, allocatable :: A(:,:), B(:,:) double precision :: tmp_re, tmp_im - integer :: n_integrals - integer(key_kind), allocatable :: buffer_i(:) - real(integral_kind), allocatable :: buffer_values(:) + integer :: n_integrals_1, n_integrals_2 + integer(key_kind), allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind), allocatable :: buffer_values_1(:), buffer_values_2(:) + logical :: use_map1 + integer(key_kind) :: idx_tmp + double precision :: sign + call ezfio_set_ao_basis_ao_num(ao_num) @@ -107,31 +111,66 @@ subroutine run call ezfio_set_ao_one_e_ints_ao_integrals_n_e_imag(B(1:ao_num, 1:ao_num)) call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e("Read") -! allocate(buffer_i(ao_num**3), buffer_values(ao_num**3)) -! iunit = getunitandopen('W.qp','r') -! n_integrals=0 -! buffer_values = 0.d0 -! do -! read (iunit,*,end=13) i,j,k,l, integral -! n_integrals += 1 -! call two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) ) -! buffer_values(n_integrals) = integral -! if (n_integrals == size(buffer_i)) then -! call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values) -! n_integrals = 0 -! endif -! enddo -! 13 continue -! close(iunit) -! -! if (n_integrals > 0) then -! call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values) -! endif -! -! call map_sort(ao_integrals_map) -! call map_unique(ao_integrals_map) -! -! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) -! call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + allocate(buffer_i_1(ao_num**3), buffer_values_1(ao_num**3)) + allocate(buffer_i_2(ao_num**3), buffer_values_2(ao_num**3)) + iunit = getunitandopen('W.qp','r') + n_integrals_1=0 + n_integrals_2=0 + buffer_values_1 = 0.d0 + buffer_values_2 = 0.d0 + do + read (iunit,*,end=13) i,j,k,l, tmp_re, tmp_im + call ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) + print'(4(I4),(L3),(I6),(F7.1))',i,j,k,l,use_map1,idx_tmp,sign + if (use_map1) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp + buffer_values_1(n_integrals_1)=tmp_re + print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp,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 + print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp+1,tmp_im*sign + endif + if (n_integrals_1 >= size(buffer_i_1)-1) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + 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 + print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp,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 + print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp+1,tmp_im*sign + endif + if (n_integrals_2 >= size(buffer_i_2)-1) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + n_integrals_2 = 0 + endif + endif + enddo + 13 continue + close(iunit) + + if (n_integrals_1 > 0) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + endif + if (n_integrals_2 > 0) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + endif + + call map_sort(ao_integrals_map) + call map_unique(ao_integrals_map) + call map_sort(ao_integrals_map_2) + call map_unique(ao_integrals_map_2) + + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_1',ao_integrals_map) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_2',ao_integrals_map_2) + call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') end