From 451351600e75d1262b69190ec8368029febc51fc Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 25 Apr 2025 11:32:33 +0200 Subject: [PATCH] problem with python --- src/ao_two_e_ints/EZFIO.cfg | 19 ++ src/ao_two_e_ints/map_integrals.irp.f | 368 +++++++++------------- src/ao_two_e_ints/map_integrals_erf.irp.f | 288 ----------------- src/ao_two_e_ints/periodic.irp.f | 137 ++++++++ src/two_e_ints_keywords/direct.irp.f | 4 + 5 files changed, 317 insertions(+), 499 deletions(-) delete mode 100644 src/ao_two_e_ints/map_integrals_erf.irp.f create mode 100644 src/ao_two_e_ints/periodic.irp.f create mode 100644 src/two_e_ints_keywords/direct.irp.f diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index b6fd2489..149e0733 100644 --- a/src/ao_two_e_ints/EZFIO.cfg +++ b/src/ao_two_e_ints/EZFIO.cfg @@ -15,3 +15,22 @@ type: Disk_access doc: Read/Write |AO| erf integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None + +[io_ao_erf_cholesky] +type: Disk_access +doc: Read/Write |AO| erf Cholesky integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[io_ao_two_e_integrals_cgtos] +type: Disk_access +doc: Read/Write |AO| cgtos integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[io_ao_cgtos_cholesky] +type: Disk_access +doc: Read/Write |AO| cgtos Cholesky integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index e429047f..17a01f44 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -2,39 +2,38 @@ use map_module !! AO Map !! ====== +BEGIN_TEMPLATE -BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] +BEGIN_PROVIDER [ type(map_type), ao_integrals$_erf_map ] implicit none BEGIN_DOC - ! AO integrals + ! AO integrals$_erf END_DOC 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) + call map_init(ao_integrals$_erf_map,sze) print*, 'AO map initialized : ', sze END_PROVIDER - - - BEGIN_PROVIDER [ integer, ao_integrals_cache_min ] -&BEGIN_PROVIDER [ integer, ao_integrals_cache_max ] + BEGIN_PROVIDER [ integer, ao_integrals$_erf_cache_min ] +&BEGIN_PROVIDER [ integer, ao_integrals$_erf_cache_max ] implicit none BEGIN_DOC - ! Min and max values of the AOs for which the integrals are in the cache + ! Min and max values of the AOs for which the integrals$_erf are in the cache END_DOC - ao_integrals_cache_min = max(1,ao_num - 63) - ao_integrals_cache_max = ao_num + ao_integrals$_erf_cache_min = max(1,ao_num - 63) + ao_integrals$_erf_cache_max = ao_num END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ] +BEGIN_PROVIDER [ double precision, ao_integrals$_erf_cache, (0:64*64*64*64) ] implicit none BEGIN_DOC - ! Cache of AO integrals for fast access + ! Cache of AO integrals$_erf for fast access END_DOC - PROVIDE ao_two_e_integrals_in_map + PROVIDE ao_two_e_integrals$_erf_in_map integer :: i,j,k,l,ii integer(key_kind) :: idx, idx2 real(integral_kind) :: integral @@ -42,19 +41,19 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ] integer(key_kind) :: idx_re,idx_im !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,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 + do l=ao_integrals$_erf_cache_min,ao_integrals$_erf_cache_max + do k=ao_integrals$_erf_cache_min,ao_integrals$_erf_cache_max + do j=ao_integrals$_erf_cache_min,ao_integrals$_erf_cache_max + do i=ao_integrals$_erf_cache_min,ao_integrals$_erf_cache_max !DIR$ FORCEINLINE call two_e_integrals_index(i,j,k,l,idx) !DIR$ FORCEINLINE - call map_get(ao_integrals_map,idx,integral) - 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) - ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) - ao_integrals_cache(ii) = integral + call map_get(ao_integrals$_erf_map,idx,integral) + ii = l-ao_integrals$_erf_cache_min + ii = ior( shiftl(ii,6), k-ao_integrals$_erf_cache_min) + ii = ior( shiftl(ii,6), j-ao_integrals$_erf_cache_min) + ii = ior( shiftl(ii,6), i-ao_integrals$_erf_cache_min) + ao_integrals$_erf_cache(ii) = integral enddo enddo enddo @@ -64,7 +63,7 @@ END_PROVIDER ! --- -double precision function get_ao_two_e_integral(i, j, k, l, map) result(result) +double precision function get_ao_two_e_integral$_erf(i, j, k, l) result(result) use map_module implicit none BEGIN_DOC @@ -74,143 +73,35 @@ double precision function get_ao_two_e_integral(i, j, k, l, map) result(result) END_DOC integer, intent(in) :: i,j,k,l integer(key_kind) :: idx - type(map_type), intent(inout) :: map integer :: ii real(integral_kind) :: tmp logical, external :: ao_two_e_integral_zero - PROVIDE ao_two_e_integrals_in_map ao_integrals_cache ao_integrals_cache_min + PROVIDE ao_two_e_integrals$_erf_in_map ao_integrals$_erf_cache ao_integrals$_erf_cache_min !DIR$ FORCEINLINE if (ao_two_e_integral_zero(i,j,k,l)) then tmp = 0.d0 else - ii = l-ao_integrals_cache_min - ii = ior(ii, k-ao_integrals_cache_min) - ii = ior(ii, j-ao_integrals_cache_min) - ii = ior(ii, i-ao_integrals_cache_min) + ii = l-ao_integrals$_erf_cache_min + ii = ior(ii, k-ao_integrals$_erf_cache_min) + ii = ior(ii, j-ao_integrals$_erf_cache_min) + ii = ior(ii, i-ao_integrals$_erf_cache_min) if (iand(ii, -64) /= 0) then !DIR$ FORCEINLINE call two_e_integrals_index(i,j,k,l,idx) !DIR$ FORCEINLINE call map_get(map,idx,tmp) else - 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) - ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) - tmp = ao_integrals_cache(ii) + ii = l-ao_integrals$_erf_cache_min + ii = ior( shiftl(ii,6), k-ao_integrals$_erf_cache_min) + ii = ior( shiftl(ii,6), j-ao_integrals$_erf_cache_min) + ii = ior( shiftl(ii,6), i-ao_integrals$_erf_cache_min) + tmp = ao_integrals$_erf_cache(ii) endif endif result = tmp end -BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] - implicit none - BEGIN_DOC - ! Cache of AO integrals for fast access - END_DOC - PROVIDE ao_two_e_integrals_in_map - integer :: i,j,k,l,ii - integer(key_kind) :: idx1, idx2 - real(integral_kind) :: tmp_re, tmp_im - integer(key_kind) :: idx_re,idx_im - complex(integral_kind) :: integral - - - !$OMP PARALLEL DO PRIVATE (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 - 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 - - 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) - ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) - ao_integrals_cache_periodic(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - -complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) 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 - integer :: ii - complex(integral_kind) :: tmp - PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min - !DIR$ FORCEINLINE - logical, external :: ao_two_e_integral_zero - if (ao_two_e_integral_zero(i,j,k,l)) then - tmp = (0.d0,0.d0) - else - ii = l-ao_integrals_cache_min - ii = ior(ii, k-ao_integrals_cache_min) - 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 - else - 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) - ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) - tmp = ao_integrals_cache_periodic(ii) - endif - result = tmp - endif -end - - -subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) +subroutine get_ao_two_e_integrals$_erf(j,k,l,sze,out_val) use map_module BEGIN_DOC ! Gets multiple AO bi-electronic integral from the AO map . @@ -224,50 +115,22 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) integer :: i integer(key_kind) :: hash logical, external :: ao_one_e_integral_zero - PROVIDE ao_two_e_integrals_in_map ao_integrals_map + PROVIDE ao_two_e_integrals$_erf_in_map ao_integrals$_erf_map if (ao_one_e_integral_zero(j,l)) then out_val(1:sze) = 0.d0 return endif - double precision :: get_ao_two_e_integral + double precision :: get_ao_two_e_integral$_erf do i=1,sze - out_val(i) = get_ao_two_e_integral(i,j,k,l,ao_integrals_map) + out_val(i) = get_ao_two_e_integral$_erf(i,j,k,l) enddo end -subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val) - use map_module - BEGIN_DOC - ! Gets multiple AO bi-electronic integral from the AO map . - ! All i are retrieved for j,k,l fixed. - ! physicist convention : - END_DOC - implicit none - integer, intent(in) :: j,k,l, sze - complex(integral_kind), intent(out) :: out_val(sze) - - integer :: i - integer(key_kind) :: hash - logical, external :: ao_one_e_integral_zero - PROVIDE ao_two_e_integrals_in_map ao_integrals_map - - if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - double precision :: get_ao_two_e_integral - do i=1,sze - out_val(i) = get_ao_two_e_integral(i,j,k,l,ao_integrals_map) - enddo - -end - -subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) +subroutine get_ao_two_e_integrals$_erf_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) use map_module implicit none BEGIN_DOC @@ -283,24 +146,23 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z double precision :: tmp logical, external :: ao_one_e_integral_zero logical, external :: ao_two_e_integral_zero - PROVIDE ao_two_e_integrals_in_map + PROVIDE ao_two_e_integrals$_erf_in_map - non_zero_int = 0 - if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif + !non_zero_int = 0 + !if (ao_one_e_integral_zero(j,l)) then + ! out_val = 0.d0 + ! return + !endif non_zero_int = 0 do i=1,sze integer, external :: ao_l4 - double precision, external :: ao_two_e_integral !DIR$ FORCEINLINE - if (ao_two_e_integral_zero(i,j,k,l)) then - cycle - endif + !if (ao_two_e_integral_zero(i,j,k,l)) then + ! cycle + !endif call two_e_integrals_index(i,j,k,l,hash) - call map_get(ao_integrals_map, hash,tmp) + call map_get(ao_integrals$_erf_map, hash,tmp) if (dabs(tmp) < ao_integrals_threshold) cycle non_zero_int = non_zero_int+1 out_val_index(non_zero_int) = i @@ -310,7 +172,7 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z end -subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int) +subroutine get_ao_two_e_integrals$_erf_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int) use map_module implicit none BEGIN_DOC @@ -328,12 +190,12 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out logical, external :: ao_one_e_integral_zero logical, external :: ao_two_e_integral_zero - PROVIDE ao_two_e_integrals_in_map - non_zero_int = 0 - if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif + PROVIDE ao_two_e_integrals$_erf_in_map + !non_zero_int = 0 + !if (ao_one_e_integral_zero(j,l)) then + ! out_val = 0.d0 + ! return + !endif non_zero_int = 0 do k = 1, sze @@ -345,7 +207,7 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out cycle endif call two_e_integrals_index(i,j,k,l,hash) - call map_get(ao_integrals_map, hash,tmp) + call map_get(ao_integrals$_erf_map, hash,tmp) if (dabs(tmp) < thresh ) cycle non_zero_int = non_zero_int+1 out_val_index(1,non_zero_int) = i @@ -357,11 +219,11 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out end -subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int) +subroutine get_ao_two_e_integrals$_erf_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int) use map_module implicit none BEGIN_DOC - ! Gets multiple AO two-electron integrals from the AO map . + ! Gets multiple AO two-electron integrals$_erf from the AO map . ! All non-zero i are retrieved for j,k,l fixed. END_DOC double precision, intent(in) :: thresh @@ -373,15 +235,15 @@ subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,s integer :: i,k integer(key_kind) :: hash double precision :: tmp - logical, external :: ao_one_e_integral_zero - logical, external :: ao_two_e_integral_zero +! logical, external :: ao_one_e_integral_zero +! logical, external :: ao_two_e_integral_zero - PROVIDE ao_two_e_integrals_in_map + PROVIDE ao_two_e_integrals$_erf_in_map non_zero_int = 0 - if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif +! if (ao_one_e_integral_zero(j,l)) then +! out_val = 0.d0 +! return +! endif non_zero_int = 0 integer :: kk @@ -395,7 +257,7 @@ subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,s cycle endif call two_e_integrals_index(i,j,k,l,hash) - call map_get(ao_integrals_map, hash,tmp) + call map_get(ao_integrals$_erf_map, hash,tmp) if (dabs(tmp) < thresh ) cycle non_zero_int = non_zero_int+1 out_val_index(1,non_zero_int) = i @@ -408,26 +270,26 @@ end -function get_ao_map_size() +function get_ao$_erf_map_size() implicit none integer (map_size_kind) :: 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$_erf_size = ao_integrals$_erf_map % n_elements end -subroutine clear_ao_map +subroutine clear_ao$_erf_map implicit none BEGIN_DOC ! Frees the memory of the AO map END_DOC - call map_deinit(ao_integrals_map) - FREE ao_integrals_map + call map_deinit(ao_integrals$_erf_map) + FREE ao_integrals$_erf_map end -subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values) +subroutine insert_into_ao_integrals$_erf_map(n_integrals,buffer_i, buffer_values) use map_module implicit none BEGIN_DOC @@ -438,7 +300,91 @@ subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values) integer(key_kind), intent(inout) :: buffer_i(n_integrals) real(integral_kind), intent(inout) :: buffer_values(n_integrals) - call map_append(ao_integrals_map, buffer_i, buffer_values, n_integrals) + call map_append(ao_integrals$_erf_map, buffer_i, buffer_values, n_integrals) +end + +subroutine dump_ao_integrals$_erf(filename) + use map_module + implicit none + BEGIN_DOC + ! Save to disk the |AO| $_erf integrals + END_DOC + character*(*), intent(in) :: filename + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer*8 :: i,j, n + call ezfio_set_work_empty(.False.) + open(unit=66,file=filename,FORM='unformatted') + write(66) integral_kind, key_kind + write(66) ao_integral$s_erf_map%sorted, ao_integrals$_erf_map%map_size, & + ao_integrals$_erf_map%n_elements + do i=0_8,ao_integrals$_erf_map%map_size + write(66) ao_integrals$_erf_map%map(i)%sorted, ao_integrals$_erf_map%map(i)%map_size,& + ao_integrals$_erf_map%map(i)%n_elements + enddo + do i=0_8,ao_integrals$_erf_map%map_size + key => ao_integrals$_erf_map%map(i)%key + val => ao_integrals$_erf_map%map(i)%value + n = ao_integrals$_erf_map%map(i)%n_elements + write(66) (key(j), j=1,n), (val(j), j=1,n) + enddo + close(66) + end + +integer function load_ao_integrals$_erf(filename) + implicit none + BEGIN_DOC + ! Read from disk the |AO| erf 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$_erf = 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$_erf_map%sorted, ao_integrals$_erf_map%map_size,& + ao_integrals$_erf_map%n_elements + do i=0_8, ao_integrals$_erf_map%map_size + read(66,err=99,end=99) ao_integrals$_erf_map%map(i)%sorted, & + ao_integrals$_erf_map%map(i)%map_size, ao_integrals$_erf_map%map(i)%n_elements + call cache_map_reallocate(ao_integrals$_erf_map%map(i),ao_integrals$_erf_map%map(i)%map_size) + enddo + do i=0_8, ao_integrals$_erf_map%map_size + key => ao_integrals$_erf_map%map(i)%key + val => ao_integrals$_erf_map%map(i)%value + n = ao_integrals$_erf_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$_erf_map) + load_ao_integrals$_erf = 0 + return + 99 continue + call map_deinit(ao_integrals$_erf_map) + 98 continue + stop 'Problem reading ao_integrals$_erf_map file in work/' + +end + + + +SUBST [ _erf ] + +;; +_erf;; +_cgtos;; + +END_TEMPLATE diff --git a/src/ao_two_e_ints/map_integrals_erf.irp.f b/src/ao_two_e_ints/map_integrals_erf.irp.f deleted file mode 100644 index 4e4e68ff..00000000 --- a/src/ao_two_e_ints/map_integrals_erf.irp.f +++ /dev/null @@ -1,288 +0,0 @@ -use map_module - -!! AO Map -!! ====== - -BEGIN_PROVIDER [ type(map_type), ao_integrals_erf_map ] - implicit none - BEGIN_DOC - ! |AO| integrals - END_DOC - 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_erf_map,sze) - print*, 'AO map initialized : ', sze -END_PROVIDER - - BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_min ] -&BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_max ] - implicit none - BEGIN_DOC - ! Min and max values of the AOs for which the integrals are in the cache - END_DOC - ao_integrals_erf_cache_min = max(1,ao_num - 63) - ao_integrals_erf_cache_max = ao_num - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ] - use map_module - implicit none - BEGIN_DOC - ! Cache of |AO| integrals for fast access - END_DOC - PROVIDE ao_two_e_integrals_erf_in_map - integer :: i,j,k,l,ii - integer(key_kind) :: idx - real(integral_kind) :: integral - !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) - do l=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max - do k=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max - do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max - do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx) - !DIR$ FORCEINLINE - call map_get(ao_integrals_erf_map,idx,integral) - ii = l-ao_integrals_erf_cache_min - ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min) - ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min) - ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min) - ao_integrals_erf_cache(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - -subroutine insert_into_ao_integrals_erf_map(n_integrals,buffer_i, buffer_values) - use map_module - implicit none - BEGIN_DOC - ! Create new entry into |AO| map - 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) - - call map_append(ao_integrals_erf_map, buffer_i, buffer_values, n_integrals) -end - -double precision function get_ao_two_e_integral_erf(i,j,k,l,map) result(result) - use map_module - implicit none - BEGIN_DOC - ! Gets one |AO| two-electron integral from the |AO| map - END_DOC - integer, intent(in) :: i,j,k,l - integer(key_kind) :: idx - type(map_type), intent(inout) :: map - integer :: ii - real(integral_kind) :: tmp - logical, external :: ao_two_e_integral_zero - PROVIDE ao_two_e_integrals_erf_in_map ao_integrals_erf_cache ao_integrals_erf_cache_min - !DIR$ FORCEINLINE - if (ao_two_e_integral_zero(i,j,k,l)) then - tmp = 0.d0 - else if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < ao_integrals_threshold) then - tmp = 0.d0 - else - ii = l-ao_integrals_erf_cache_min - ii = ior(ii, k-ao_integrals_erf_cache_min) - ii = ior(ii, j-ao_integrals_erf_cache_min) - ii = ior(ii, i-ao_integrals_erf_cache_min) - if (iand(ii, -64) /= 0) then - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx) - !DIR$ FORCEINLINE - call map_get(map,idx,tmp) - tmp = tmp - else - ii = l-ao_integrals_erf_cache_min - ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min) - ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min) - ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min) - tmp = ao_integrals_erf_cache(ii) - endif - endif - result = tmp -end - - -subroutine get_ao_two_e_integrals_erf(j,k,l,sze,out_val) - use map_module - BEGIN_DOC - ! Gets multiple |AO| two-electron integral from the |AO| map . - ! All i are retrieved for j,k,l fixed. - END_DOC - implicit none - integer, intent(in) :: j,k,l, sze - real(integral_kind), intent(out) :: out_val(sze) - - integer :: i - integer(key_kind) :: hash - double precision :: thresh - logical, external :: ao_one_e_integral_zero - PROVIDE ao_two_e_integrals_erf_in_map ao_integrals_erf_map - thresh = ao_integrals_threshold - - if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - double precision :: get_ao_two_e_integral_erf - do i=1,sze - out_val(i) = get_ao_two_e_integral_erf(i,j,k,l,ao_integrals_erf_map) - enddo - -end - -subroutine get_ao_two_e_integrals_erf_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) - use map_module - implicit none - BEGIN_DOC - ! Gets multiple |AO| two-electron integrals from the |AO| map . - ! All non-zero i are retrieved for j,k,l fixed. - END_DOC - integer, intent(in) :: j,k,l, sze - real(integral_kind), intent(out) :: out_val(sze) - integer, intent(out) :: out_val_index(sze),non_zero_int - - integer :: i - integer(key_kind) :: hash - double precision :: thresh,tmp - logical, external :: ao_one_e_integral_zero - PROVIDE ao_two_e_integrals_erf_in_map - thresh = ao_integrals_threshold - - non_zero_int = 0 - if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - non_zero_int = 0 - do i=1,sze - integer, external :: ao_l4 - double precision, external :: ao_two_e_integral_erf - !DIR$ FORCEINLINE - if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thresh) then - cycle - endif - call two_e_integrals_index(i,j,k,l,hash) - call map_get(ao_integrals_erf_map, hash,tmp) - if (dabs(tmp) < thresh ) cycle - non_zero_int = non_zero_int+1 - out_val_index(non_zero_int) = i - out_val(non_zero_int) = tmp - enddo - -end - - -function get_ao_erf_map_size() - implicit none - integer (map_size_kind) :: get_ao_erf_map_size - BEGIN_DOC - ! Returns the number of elements in the |AO| map - END_DOC - get_ao_erf_map_size = ao_integrals_erf_map % n_elements -end - -subroutine clear_ao_erf_map - implicit none - BEGIN_DOC - ! Frees the memory of the |AO| map - END_DOC - call map_deinit(ao_integrals_erf_map) - FREE ao_integrals_erf_map -end - - - -subroutine dump_ao_integrals_erf(filename) - use map_module - implicit none - BEGIN_DOC - ! Save to disk the |AO| erf integrals - END_DOC - character*(*), intent(in) :: filename - integer(cache_key_kind), pointer :: key(:) - real(integral_kind), pointer :: val(:) - integer*8 :: i,j, n - call ezfio_set_work_empty(.False.) - open(unit=66,file=filename,FORM='unformatted') - write(66) integral_kind, key_kind - write(66) ao_integrals_erf_map%sorted, ao_integrals_erf_map%map_size, & - ao_integrals_erf_map%n_elements - do i=0_8,ao_integrals_erf_map%map_size - write(66) ao_integrals_erf_map%map(i)%sorted, ao_integrals_erf_map%map(i)%map_size,& - ao_integrals_erf_map%map(i)%n_elements - enddo - do i=0_8,ao_integrals_erf_map%map_size - key => ao_integrals_erf_map%map(i)%key - val => ao_integrals_erf_map%map(i)%value - n = ao_integrals_erf_map%map(i)%n_elements - write(66) (key(j), j=1,n), (val(j), j=1,n) - enddo - close(66) - -end - - - -integer function load_ao_integrals_erf(filename) - implicit none - BEGIN_DOC - ! Read from disk the |AO| erf 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_erf = 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_erf_map%sorted, ao_integrals_erf_map%map_size,& - ao_integrals_erf_map%n_elements - do i=0_8, ao_integrals_erf_map%map_size - read(66,err=99,end=99) ao_integrals_erf_map%map(i)%sorted, & - ao_integrals_erf_map%map(i)%map_size, ao_integrals_erf_map%map(i)%n_elements - call cache_map_reallocate(ao_integrals_erf_map%map(i),ao_integrals_erf_map%map(i)%map_size) - enddo - do i=0_8, ao_integrals_erf_map%map_size - key => ao_integrals_erf_map%map(i)%key - val => ao_integrals_erf_map%map(i)%value - n = ao_integrals_erf_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_erf_map) - load_ao_integrals_erf = 0 - return - 99 continue - call map_deinit(ao_integrals_erf_map) - 98 continue - stop 'Problem reading ao_integrals_erf_map file in work/' - -end - - - - diff --git a/src/ao_two_e_ints/periodic.irp.f b/src/ao_two_e_ints/periodic.irp.f new file mode 100644 index 00000000..a13c04a5 --- /dev/null +++ b/src/ao_two_e_ints/periodic.irp.f @@ -0,0 +1,137 @@ + +BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] + implicit none + BEGIN_DOC + ! Cache of AO integrals for fast access + END_DOC + PROVIDE ao_two_e_integrals_in_map + integer :: i,j,k,l,ii + integer(key_kind) :: idx1, idx2 + real(integral_kind) :: tmp_re, tmp_im + integer(key_kind) :: idx_re,idx_im + complex(integral_kind) :: integral + + + !$OMP PARALLEL DO PRIVATE (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 + 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 + + 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) + ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) + ao_integrals_cache_periodic(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) 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 + integer :: ii + complex(integral_kind) :: tmp + PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min + !DIR$ FORCEINLINE + logical, external :: ao_two_e_integral_zero + if (ao_two_e_integral_zero(i,j,k,l)) then + tmp = (0.d0,0.d0) + else + ii = l-ao_integrals_cache_min + ii = ior(ii, k-ao_integrals_cache_min) + 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 + else + 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) + ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) + tmp = ao_integrals_cache_periodic(ii) + endif + result = tmp + endif +end + + + +subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val) + use map_module + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All i are retrieved for j,k,l fixed. + ! physicist convention : + END_DOC + implicit none + integer, intent(in) :: j,k,l, sze + complex(integral_kind), intent(out) :: out_val(sze) + + integer :: i + integer(key_kind) :: hash + logical, external :: ao_one_e_integral_zero + PROVIDE ao_two_e_integrals_in_map ao_integrals_map + + if (ao_one_e_integral_zero(j,l)) then + out_val = 0.d0 + return + endif + + double precision :: get_ao_two_e_integral + do i=1,sze + out_val(i) = get_ao_two_e_integral(i,j,k,l,ao_integrals_map) + enddo + +end + diff --git a/src/two_e_ints_keywords/direct.irp.f b/src/two_e_ints_keywords/direct.irp.f new file mode 100644 index 00000000..9eef711a --- /dev/null +++ b/src/two_e_ints_keywords/direct.irp.f @@ -0,0 +1,4 @@ +BEGIN_PROVIDER [ logical, do_direct_integrals] + implicit none + do_direct_integrals = do_ao_cholesky +END_PROVIDER