mirror of
https://github.com/QuantumPackage/qp2.git
synced 2025-04-25 17:54:44 +02:00
big work
This commit is contained in:
eginer
parent
afe0baa60a
commit
0c76c27440
@ -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
|
||||
|
||||
@ -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 (i<j) then
|
||||
ij=(j-1)*(j-1)+2*i-mod(j+1,2)
|
||||
else if (i>j) 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 : <ij|kl>
|
||||
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 : <ij|kl>
|
||||
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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -3,3 +3,5 @@ ao_one_e_ints
|
||||
pseudo
|
||||
bitmask
|
||||
ao_basis
|
||||
two_e_ints_keywords
|
||||
ao_cart_two_e_ints
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user