From 0c76c2744037c7b36102ee7a100c43bf82929a6a Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 18 Apr 2025 16:37:44 +0200 Subject: [PATCH] big work --- src/ao_cart_two_e_ints/cholesky.irp.f | 68 +- src/ao_cart_two_e_ints/map_integrals.irp.f | 703 ------------------- src/ao_cart_two_e_ints/two_e_integrals.irp.f | 138 ---- src/ao_two_e_ints/NEED | 2 + src/ao_two_e_ints/direct.irp.f | 11 - src/ao_two_e_ints/gauss_legendre.irp.f | 57 -- src/two_e_ints_keywords/direct.irp.f | 11 - 7 files changed, 22 insertions(+), 968 deletions(-) delete mode 100644 src/ao_cart_two_e_ints/map_integrals.irp.f delete mode 100644 src/ao_two_e_ints/direct.irp.f delete mode 100644 src/ao_two_e_ints/gauss_legendre.irp.f delete mode 100644 src/two_e_ints_keywords/direct.irp.f diff --git a/src/ao_cart_two_e_ints/cholesky.irp.f b/src/ao_cart_two_e_ints/cholesky.irp.f index 8de035b6..2f68d11a 100644 --- a/src/ao_cart_two_e_ints/cholesky.irp.f +++ b/src/ao_cart_two_e_ints/cholesky.irp.f @@ -100,13 +100,9 @@ END_PROVIDER call set_multiple_levels_omp(.False.) - if (do_direct_integrals) then - if (ao_cart_two_e_integral(1,1,1,1) < huge(1.d0)) then - ! Trigger providers inside ao_cart_two_e_integral - continue - endif - else - PROVIDE ao_cart_two_e_integrals_in_map + if (ao_cart_two_e_integral(1,1,1,1) < huge(1.d0)) then + ! Trigger providers inside ao_cart_two_e_integral + continue endif tau = ao_cart_cholesky_threshold @@ -140,21 +136,12 @@ END_PROVIDER enddo enddo - if (do_direct_integrals) then - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) - do i8=ndim8-1,1,-1 - D(i8) = ao_cart_two_e_integral(addr1(i8), addr2(i8), & - addr1(i8), addr2(i8)) - enddo - !$OMP END PARALLEL DO - else - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) - do i8=ndim8-1,1,-1 - D(i8) = get_ao_cart_two_e_integral(addr1(i8), addr1(i8), & - addr2(i8), addr2(i8), ao_cart_integrals_map) - enddo - !$OMP END PARALLEL DO - endif + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) + do i8=ndim8-1,1,-1 + D(i8) = ao_cart_two_e_integral(addr1(i8), addr2(i8), & + addr1(i8), addr2(i8)) + enddo + !$OMP END PARALLEL DO ! Just to guarentee termination D(ndim8) = 0.d0 @@ -352,32 +339,17 @@ END_PROVIDER if (.not.computed(dj)) then m = dj - if (do_direct_integrals) then - !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21) - do k=1,np - Delta_col(k) = 0.d0 - if (.not.ao_cart_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& - addr2(Lset(k)), addr2(Dset(m)) ) ) then - Delta_col(k) = & - ao_cart_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& - addr1(Dset(m)), addr2(Dset(m))) - endif - enddo - !$OMP END PARALLEL DO - else - PROVIDE ao_cart_integrals_map - !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21) - do k=1,np - Delta_col(k) = 0.d0 - if (.not.ao_cart_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& - addr2(Lset(k)), addr2(Dset(m)) ) ) then - Delta_col(k) = & - get_ao_cart_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),& - addr2(Lset(k)), addr2(Dset(m)), ao_cart_integrals_map) - endif - enddo - !$OMP END PARALLEL DO - endif + !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21) + do k=1,np + Delta_col(k) = 0.d0 + if (.not.ao_cart_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& + addr2(Lset(k)), addr2(Dset(m)) ) ) then + Delta_col(k) = & + ao_cart_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& + addr1(Dset(m)), addr2(Dset(m))) + endif + enddo + !$OMP END PARALLEL DO !$OMP PARALLEL DO PRIVATE(p) do p=1,np diff --git a/src/ao_cart_two_e_ints/map_integrals.irp.f b/src/ao_cart_two_e_ints/map_integrals.irp.f deleted file mode 100644 index bf0791f1..00000000 --- a/src/ao_cart_two_e_ints/map_integrals.irp.f +++ /dev/null @@ -1,703 +0,0 @@ -use map_module - -!! AO Map -!! ====== - -BEGIN_PROVIDER [ type(map_type), ao_cart_integrals_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_cart_num,ao_cart_num,ao_cart_num,ao_cart_num,key_max) - sze = key_max - call map_init(ao_cart_integrals_map,sze) - print*, 'AO map initialized : ', sze -END_PROVIDER - -subroutine two_e_integrals_index(i,j,k,l,i1) - 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) for non-periodic systems - END_DOC - integer, intent(in) :: i,j,k,l - integer(key_kind), intent(out) :: i1 - integer(key_kind) :: p,q,r,s,i2 - 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) - use map_module - implicit none - BEGIN_DOC -! Computes the 4 indices $i,j,k,l$ from a unique index $i_1$. -! For 2 indices $i,j$ and $i \le j$, we have -! $p = i(i-1)/2 + j$. -! The key point is that because $j < i$, -! $i(i-1)/2 < p \le i(i+1)/2$. So $i$ can be found by solving -! $i^2 - i - 2p=0$. One obtains $i=1 + \sqrt{1+8p}/2$ -! and $j = p - i(i-1)/2$. -! This rule is applied 3 times. First for the symmetry of the -! pairs (i,k) and (j,l), and then for the symmetry within each pair. - END_DOC - integer, intent(out) :: i(8),j(8),k(8),l(8) - integer(key_kind), intent(in) :: i1 - integer(key_kind) :: i2,i3 - i = 0 - i2 = ceiling(0.5d0*(dsqrt(dble(shiftl(i1,3)+1))-1.d0)) - l(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i2,3)+1))-1.d0)) - i3 = i1 - shiftr(i2*i2-i2,1) - k(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i3,3)+1))-1.d0)) - j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4) - i(1) = int(i3 - shiftr(k(1)*k(1)-k(1),1),4) - - !ijkl - i(2) = i(1) !ilkj - j(2) = l(1) - k(2) = k(1) - l(2) = j(1) - - i(3) = k(1) !kjil - j(3) = j(1) - k(3) = i(1) - l(3) = l(1) - - i(4) = k(1) !klij - j(4) = l(1) - k(4) = i(1) - l(4) = j(1) - - i(5) = j(1) !jilk - j(5) = i(1) - k(5) = l(1) - l(5) = k(1) - - i(6) = j(1) !jkli - j(6) = k(1) - k(6) = l(1) - l(6) = i(1) - - i(7) = l(1) !lijk - j(7) = i(1) - k(7) = j(1) - l(7) = k(1) - - i(8) = l(1) !lkji - j(8) = k(1) - k(8) = j(1) - l(8) = i(1) - - integer :: ii, jj - do ii=2,8 - do jj=1,ii-1 - if ( (i(ii) == i(jj)).and. & - (j(ii) == j(jj)).and. & - (k(ii) == k(jj)).and. & - (l(ii) == l(jj)) ) then - i(ii) = 0 - exit - endif - enddo - enddo -! This has been tested with up to 1000 AOs, and all the reverse indices are -! correct ! We can remove the test -! do ii=1,8 -! if (i(ii) /= 0) then -! call two_e_integrals_index(i(ii),j(ii),k(ii),l(ii),i2) -! if (i1 /= i2) then -! print *, i1, i2 -! print *, i(ii), j(ii), k(ii), l(ii) -! stop 'two_e_integrals_index_reverse failed' -! endif -! endif -! enddo - - -end - - - -subroutine ao_cart_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_cart_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_cart_idx2_sq(i,k,ik) - call ao_cart_idx2_sq(j,l,jl) - call ao_cart_idx2_tri_key(ik,jl,i1) -end - -subroutine ao_cart_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_cart_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_cart_idx2_tri_rev_key(i3,i2,i1) - - call ao_cart_idx2_sq_rev(j(1),l(1),i2) - call ao_cart_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_cart_integrals_cache_min ] -&BEGIN_PROVIDER [ integer, ao_cart_integrals_cache_max ] - implicit none - BEGIN_DOC - ! Min and max values of the AOs for which the integrals are in the cache - END_DOC - ao_cart_integrals_cache_min = max(1,ao_cart_num - 63) - ao_cart_integrals_cache_max = ao_cart_num - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, ao_cart_integrals_cache, (0:64*64*64*64) ] - implicit none - BEGIN_DOC - ! Cache of AO integrals for fast access - END_DOC - PROVIDE ao_cart_two_e_integrals_in_map - integer :: i,j,k,l,ii - integer(key_kind) :: idx, idx2 - real(integral_kind) :: integral - real(integral_kind) :: tmp_re, tmp_im - integer(key_kind) :: idx_re,idx_im - - !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) - do l=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max - do k=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max - do j=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max - do i=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx) - !DIR$ FORCEINLINE - call map_get(ao_cart_integrals_map,idx,integral) - ii = l-ao_cart_integrals_cache_min - ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min) - ao_cart_integrals_cache(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO -END_PROVIDER - -! --- - -double precision function get_ao_cart_two_e_integral(i, j, k, l, map) result(result) - use map_module - implicit none - BEGIN_DOC - ! Gets one AO bi-electronic integral from the AO map in PHYSICIST NOTATION - ! - ! <1:k, 2:l |1:i, 2:j> - 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_cart_two_e_integral_zero - PROVIDE ao_cart_two_e_integrals_in_map ao_cart_integrals_cache ao_cart_integrals_cache_min - !DIR$ FORCEINLINE - if (ao_cart_two_e_integral_zero(i,j,k,l)) then - tmp = 0.d0 - else - ii = l-ao_cart_integrals_cache_min - ii = ior(ii, k-ao_cart_integrals_cache_min) - ii = ior(ii, j-ao_cart_integrals_cache_min) - ii = ior(ii, i-ao_cart_integrals_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_cart_integrals_cache_min - ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min) - tmp = ao_cart_integrals_cache(ii) - endif - endif - result = tmp -end - -BEGIN_PROVIDER [ complex*16, ao_cart_integrals_cache_periodic, (0:64*64*64*64) ] - implicit none - BEGIN_DOC - ! Cache of AO integrals for fast access - END_DOC - PROVIDE ao_cart_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_cart_integrals_cache_min,ao_cart_integrals_cache_max - do k=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max - do j=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max - do i=ao_cart_integrals_cache_min,ao_cart_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_cart_integrals_map,idx_re,tmp_re) - if (idx_re /= idx_im) then - call map_get(ao_cart_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_cart_integrals_cache_min - ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min) - ao_cart_integrals_cache_periodic(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - -complex*16 function get_ao_cart_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_cart_two_e_integrals_in_map ao_cart_integrals_cache_periodic ao_cart_integrals_cache_min - !DIR$ FORCEINLINE - logical, external :: ao_cart_two_e_integral_zero - if (ao_cart_two_e_integral_zero(i,j,k,l)) then - tmp = (0.d0,0.d0) - else - ii = l-ao_cart_integrals_cache_min - ii = ior(ii, k-ao_cart_integrals_cache_min) - ii = ior(ii, j-ao_cart_integrals_cache_min) - ii = ior(ii, i-ao_cart_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_cart_integrals_map,idx_re,tmp_re) - if (idx_re /= idx_im) then - call map_get(ao_cart_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_cart_integrals_cache_min - ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min) - ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min) - tmp = ao_cart_integrals_cache_periodic(ii) - endif - result = tmp - endif -end - - -subroutine get_ao_cart_two_e_integrals(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 - real(integral_kind), intent(out) :: out_val(sze) - - integer :: i - integer(key_kind) :: hash - logical, external :: ao_cart_one_e_integral_zero - PROVIDE ao_cart_two_e_integrals_in_map ao_cart_integrals_map - - if (ao_cart_one_e_integral_zero(j,l)) then - out_val(1:sze) = 0.d0 - return - endif - - double precision :: get_ao_cart_two_e_integral - do i=1,sze - out_val(i) = get_ao_cart_two_e_integral(i,j,k,l,ao_cart_integrals_map) - enddo - -end - - -subroutine get_ao_cart_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_cart_one_e_integral_zero - PROVIDE ao_cart_two_e_integrals_in_map ao_cart_integrals_map - - if (ao_cart_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - double precision :: get_ao_cart_two_e_integral - do i=1,sze - out_val(i) = get_ao_cart_two_e_integral(i,j,k,l,ao_cart_integrals_map) - enddo - -end - -subroutine get_ao_cart_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) - use map_module - implicit none - BEGIN_DOC - ! Gets multiple AO bi-electronic integral 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 :: tmp - logical, external :: ao_cart_one_e_integral_zero - logical, external :: ao_cart_two_e_integral_zero - PROVIDE ao_cart_two_e_integrals_in_map - - non_zero_int = 0 - if (ao_cart_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - non_zero_int = 0 - do i=1,sze - integer, external :: ao_cart_l4 - double precision, external :: ao_cart_two_e_integral - !DIR$ FORCEINLINE - if (ao_cart_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_cart_integrals_map, hash,tmp) - if (dabs(tmp) < ao_cart_integrals_threshold) cycle - non_zero_int = non_zero_int+1 - out_val_index(non_zero_int) = i - out_val(non_zero_int) = tmp - enddo - -end - - -subroutine get_ao_cart_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int) - use map_module - implicit none - BEGIN_DOC - ! Gets multiple AO bi-electronic integral from the AO map . - ! All non-zero i are retrieved for j,k,l fixed. - END_DOC - double precision, intent(in) :: thresh - integer, intent(in) :: j,l, sze,sze_max - real(integral_kind), intent(out) :: out_val(sze_max) - integer, intent(out) :: out_val_index(2,sze_max),non_zero_int - - integer :: i,k - integer(key_kind) :: hash - double precision :: tmp - logical, external :: ao_cart_one_e_integral_zero - logical, external :: ao_cart_two_e_integral_zero - - PROVIDE ao_cart_two_e_integrals_in_map - non_zero_int = 0 - if (ao_cart_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - non_zero_int = 0 - do k = 1, sze - do i = 1, sze - integer, external :: ao_cart_l4 - double precision, external :: ao_cart_two_e_integral - !DIR$ FORCEINLINE - if (ao_cart_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_cart_integrals_map, hash,tmp) - if (dabs(tmp) < thresh ) cycle - non_zero_int = non_zero_int+1 - out_val_index(1,non_zero_int) = i - out_val_index(2,non_zero_int) = k - out_val(non_zero_int) = tmp - enddo - enddo - -end - - -subroutine get_ao_cart_two_e_integrals_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 . - ! All non-zero i are retrieved for j,k,l fixed. - END_DOC - double precision, intent(in) :: thresh - integer, intent(in) :: sze_max - integer, intent(in) :: j,l, n_list,list(2,sze_max) - real(integral_kind), intent(out) :: out_val(sze_max) - integer, intent(out) :: out_val_index(2,sze_max),non_zero_int - - integer :: i,k - integer(key_kind) :: hash - double precision :: tmp - logical, external :: ao_cart_one_e_integral_zero - logical, external :: ao_cart_two_e_integral_zero - - PROVIDE ao_cart_two_e_integrals_in_map - non_zero_int = 0 - if (ao_cart_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - non_zero_int = 0 - integer :: kk - do kk = 1, n_list - k = list(1,kk) - i = list(2,kk) - integer, external :: ao_cart_l4 - double precision, external :: ao_cart_two_e_integral - !DIR$ FORCEINLINE - if (ao_cart_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_cart_integrals_map, hash,tmp) - if (dabs(tmp) < thresh ) cycle - non_zero_int = non_zero_int+1 - out_val_index(1,non_zero_int) = i - out_val_index(2,non_zero_int) = k - out_val(non_zero_int) = tmp - enddo - -end - - - - -function get_ao_cart_map_size() - implicit none - integer (map_size_kind) :: get_ao_cart_map_size - BEGIN_DOC - ! Returns the number of elements in the AO map - END_DOC - get_ao_cart_map_size = ao_cart_integrals_map % n_elements -end - -subroutine clear_ao_cart_map - implicit none - BEGIN_DOC - ! Frees the memory of the AO map - END_DOC - call map_deinit(ao_cart_integrals_map) - FREE ao_cart_integrals_map -end - - -subroutine insert_into_ao_cart_integrals_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_cart_integrals_map, buffer_i, buffer_values, n_integrals) -end - - diff --git a/src/ao_cart_two_e_ints/two_e_integrals.irp.f b/src/ao_cart_two_e_ints/two_e_integrals.irp.f index 3df0d015..7afb3095 100644 --- a/src/ao_cart_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_cart_two_e_ints/two_e_integrals.irp.f @@ -364,7 +364,6 @@ end subroutine compute_ao_cart_two_e_integrals(j,k,l,sze,buffer_value) implicit none - use map_module BEGIN_DOC ! Compute AO 1/r12 integrals for all i and fixed j,k,l @@ -396,90 +395,6 @@ subroutine compute_ao_cart_two_e_integrals(j,k,l,sze,buffer_value) end -BEGIN_PROVIDER [ logical, ao_cart_two_e_integrals_in_map ] - implicit none - use map_module - BEGIN_DOC - ! Map of Atomic integrals - ! i(r1) j(r2) 1/r12 k(r1) l(r2) - END_DOC - - integer :: i,j,k,l - double precision :: ao_cart_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - include 'utils/constants.include.F' - - ! For integrals file - integer(key_kind),allocatable :: buffer_i(:) - integer :: size_buffer - real(integral_kind),allocatable :: buffer_value(:) - - integer :: n_integrals, rc - integer :: kk, m, j1, i1, lmax - character*(64) :: fmt - - double precision :: map_mb - PROVIDE read_ao_cart_two_e_integrals io_ao_cart_two_e_integrals ao_cart_integrals_map - - if (read_ao_cart_two_e_integrals) then - print*,'Reading the AO integrals' - call map_load_from_disk(trim(ezfio_filename)//'/work/ao_cart_ints',ao_cart_integrals_map) - print*, 'AO integrals provided' - ao_cart_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) - - if (.True.) then - ! Avoid openMP - integral = ao_cart_two_e_integral(1,1,1,1) - endif - - size_buffer = ao_cart_num*ao_cart_num - !$OMP PARALLEL DEFAULT(shared) private(j,l) & - !$OMP PRIVATE(buffer_i, buffer_value, n_integrals) - allocate(buffer_i(size_buffer), buffer_value(size_buffer)) - n_integrals = 0 - !$OMP DO COLLAPSE(1) SCHEDULE(dynamic) - do l=1,ao_cart_num - do j=1,l - call compute_ao_cart_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) - call insert_into_ao_cart_integrals_map(n_integrals,buffer_i,buffer_value) - enddo - enddo - !$OMP END DO - deallocate(buffer_i, buffer_value) - !$OMP END PARALLEL - - print*, 'Sorting the map' - call map_sort(ao_cart_integrals_map) - call cpu_time(cpu_2) - call wall_time(wall_2) - integer(map_size_kind) :: get_ao_cart_map_size, ao_cart_map_size - ao_cart_map_size = get_ao_cart_map_size() - - print*, 'AO integrals provided:' - print*, ' Size of AO map : ', map_mb(ao_cart_integrals_map) ,'MB' - print*, ' Number of AO integrals :', ao_cart_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_cart_two_e_integrals_in_map = .True. - - if (write_ao_cart_two_e_integrals.and.mpi_master) then - call ezfio_set_work_empty(.False.) - call map_save_to_disk(trim(ezfio_filename)//'/work/ao_cart_ints',ao_cart_integrals_map) - call ezfio_set_ao_cart_two_e_ints_io_ao_cart_two_e_integrals('Read') - endif - -END_PROVIDER - -! --- - BEGIN_PROVIDER [ double precision, ao_cart_two_e_integral_schwartz, (ao_cart_num, ao_cart_num) ] BEGIN_DOC @@ -1653,59 +1568,6 @@ end -subroutine compute_ao_cart_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) - implicit none - use map_module - BEGIN_DOC - ! Parallel client for AO integrals - END_DOC - - integer, intent(in) :: j,l - integer,intent(out) :: n_integrals - integer(key_kind),intent(out) :: buffer_i(ao_cart_num*ao_cart_num) - real(integral_kind),intent(out) :: buffer_value(ao_cart_num*ao_cart_num) - logical, external :: ao_cart_two_e_integral_zero - - integer :: i,k - double precision, external :: ao_cart_two_e_integral - double precision :: cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - double precision :: thr - integer :: kk, m, j1, i1 - - thr = ao_cart_integrals_threshold - - n_integrals = 0 - - j1 = j+shiftr(l*l-l,1) - do k = 1, ao_cart_num ! r1 - i1 = shiftr(k*k-k,1) - if (i1 > j1) then - exit - endif - do i = 1, k - i1 += 1 - if (i1 > j1) then - exit - endif - if (ao_cart_two_e_integral_zero(i,j,k,l)) then - cycle - endif - !DIR$ FORCEINLINE - integral = ao_cart_two_e_integral(i,k,j,l) ! i,k : r1 j,l : r2 - if (abs(integral) < thr) then - cycle - endif - n_integrals += 1 - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) - buffer_value(n_integrals) = integral - enddo - enddo - -end - - subroutine multiply_poly_local(b,nb,c,nc,d,nd) implicit none BEGIN_DOC diff --git a/src/ao_two_e_ints/NEED b/src/ao_two_e_ints/NEED index 34b4a641..622e3952 100644 --- a/src/ao_two_e_ints/NEED +++ b/src/ao_two_e_ints/NEED @@ -3,3 +3,5 @@ ao_one_e_ints pseudo bitmask ao_basis +two_e_ints_keywords +ao_cart_two_e_ints diff --git a/src/ao_two_e_ints/direct.irp.f b/src/ao_two_e_ints/direct.irp.f deleted file mode 100644 index ce093a57..00000000 --- a/src/ao_two_e_ints/direct.irp.f +++ /dev/null @@ -1,11 +0,0 @@ -BEGIN_PROVIDER [ logical, do_direct_integrals ] - implicit none - BEGIN_DOC -! Compute integrals on the fly - END_DOC - - do_direct_integrals = do_ao_cholesky - -END_PROVIDER - - diff --git a/src/ao_two_e_ints/gauss_legendre.irp.f b/src/ao_two_e_ints/gauss_legendre.irp.f deleted file mode 100644 index 4bdadb6e..00000000 --- a/src/ao_two_e_ints/gauss_legendre.irp.f +++ /dev/null @@ -1,57 +0,0 @@ - BEGIN_PROVIDER [ double precision, gauleg_t2, (n_pt_max_integrals,n_pt_max_integrals/2) ] -&BEGIN_PROVIDER [ double precision, gauleg_w, (n_pt_max_integrals,n_pt_max_integrals/2) ] - implicit none - BEGIN_DOC - ! t_w(i,1,k) = w(i) - ! t_w(i,2,k) = t(i) - END_DOC - integer :: i,j,l - l=0 - do i = 2,n_pt_max_integrals,2 - l = l+1 - call gauleg(0.d0,1.d0,gauleg_t2(1,l),gauleg_w(1,l),i) - do j=1,i - gauleg_t2(j,l) *= gauleg_t2(j,l) - enddo - enddo - -END_PROVIDER - -subroutine gauleg(x1,x2,x,w,n) - implicit none - BEGIN_DOC - ! Gauss-Legendre - END_DOC - integer, intent(in) :: n - double precision, intent(in) :: x1, x2 - double precision, intent (out) :: x(n),w(n) - double precision, parameter :: eps=3.d-14 - - integer :: m,i,j - double precision :: xm, xl, z, z1, p1, p2, p3, pp, dn - m=(n+1)/2 - xm=0.5d0*(x2+x1) - xl=0.5d0*(x2-x1) - dn = dble(n) - do i=1,m - z=dcos(3.141592654d0*(dble(i)-.25d0)/(dble(n)+.5d0)) - z1 = z+1.d0 - do while (dabs(z-z1) > eps) - p1=1.d0 - p2=0.d0 - do j=1,n - p3=p2 - p2=p1 - p1=(dble(j+j-1)*z*p2-dble(j-1)*p3)/j - enddo - pp=dn*(z*p1-p2)/(z*z-1.d0) - z1=z - z=z1-p1/pp - end do - x(i)=xm-xl*z - x(n+1-i)=xm+xl*z - w(i)=(xl+xl)/((1.d0-z*z)*pp*pp) - w(n+1-i)=w(i) - enddo -end - diff --git a/src/two_e_ints_keywords/direct.irp.f b/src/two_e_ints_keywords/direct.irp.f deleted file mode 100644 index ce093a57..00000000 --- a/src/two_e_ints_keywords/direct.irp.f +++ /dev/null @@ -1,11 +0,0 @@ -BEGIN_PROVIDER [ logical, do_direct_integrals ] - implicit none - BEGIN_DOC -! Compute integrals on the fly - END_DOC - - do_direct_integrals = do_ao_cholesky - -END_PROVIDER - -