10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-25 22:52:15 +02:00

Accelerated integral access

This commit is contained in:
Anthony Scemama 2016-10-28 21:59:39 +02:00
parent 08ac74cc2d
commit 156a3f551b

View File

@ -120,15 +120,16 @@ end
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrals_cache, (ao_integrals_cache_min:ao_integrals_cache_max,ao_integrals_cache_min:ao_integrals_cache_max,ao_integrals_cache_min:ao_integrals_cache_max,ao_integrals_cache_min:ao_integrals_cache_max) ]
BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_bielec_integrals_in_map
integer :: i,j,k,l
integer :: i,j,k,l,ii
integer(key_kind) :: idx
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx)
real(integral_kind) :: integral
!$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
@ -136,7 +137,12 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (ao_integrals_cache_min:a
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(ao_integrals_map,idx,ao_integrals_cache(i,j,k,l))
call map_get(ao_integrals_map,idx,integral)
ii = l-ao_integrals_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_cache_min)
ao_integrals_cache(ii) = integral
enddo
enddo
enddo
@ -155,30 +161,33 @@ double precision function get_ao_bielec_integral(i,j,k,l,map)
integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx
type(map_type), intent(inout) :: map
integer :: ii
real(integral_kind) :: tmp
PROVIDE ao_bielec_integrals_in_map
PROVIDE ao_bielec_integrals_in_map ao_integrals_cache ao_integrals_cache_min
!DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
tmp = 0.d0
else if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < ao_integrals_threshold) then
tmp = 0.d0
else
if ( (i >= ao_integrals_cache_min) .and. &
(j >= ao_integrals_cache_min) .and. &
(k >= ao_integrals_cache_min) .and. &
(l >= ao_integrals_cache_min) .and. &
(i <= ao_integrals_cache_max) .and. &
(j <= ao_integrals_cache_max) .and. &
(k <= ao_integrals_cache_max) .and. &
(l <= ao_integrals_cache_max) ) then
tmp = ao_integrals_cache(i,j,k,l)
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 bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
get_ao_bielec_integral = dble(tmp)
else
ii = l-ao_integrals_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_cache_min)
get_ao_bielec_integral = ao_integrals_cache(ii)
endif
endif
get_ao_bielec_integral = tmp
end
@ -324,20 +333,22 @@ end
! Min and max values of the MOs for which the integrals are in the cache
END_DOC
mo_integrals_cache_min = max(1,elec_alpha_num - 31)
mo_integrals_cache_max = min(mo_tot_num,elec_alpha_num + 32)
mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+63)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_integrals_cache, (mo_integrals_cache_min:mo_integrals_cache_max,mo_integrals_cache_min:mo_integrals_cache_max,mo_integrals_cache_min:mo_integrals_cache_max,mo_integrals_cache_min:mo_integrals_cache_max) ]
BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:64*64*64*64) ]
implicit none
BEGIN_DOC
! Cache of MO integrals for fast access
END_DOC
PROVIDE mo_bielec_integrals_in_map
integer :: i,j,k,l
integer :: ii
integer(key_kind) :: idx
real(integral_kind) :: integral
FREE ao_integrals_cache
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx)
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=mo_integrals_cache_min,mo_integrals_cache_max
do k=mo_integrals_cache_min,mo_integrals_cache_max
do j=mo_integrals_cache_min,mo_integrals_cache_max
@ -345,7 +356,12 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (mo_integrals_cache_min:m
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(mo_integrals_map,idx,mo_integrals_cache(i,j,k,l))
call map_get(mo_integrals_map,idx,integral)
ii = l-mo_integrals_cache_min
ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
mo_integrals_cache(ii) = integral
enddo
enddo
enddo
@ -362,24 +378,26 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
END_DOC
integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx
integer :: ii
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
if ( (i >= mo_integrals_cache_min) .and. &
(j >= mo_integrals_cache_min) .and. &
(k >= mo_integrals_cache_min) .and. &
(l >= mo_integrals_cache_min) .and. &
(i <= mo_integrals_cache_max) .and. &
(j <= mo_integrals_cache_max) .and. &
(k <= mo_integrals_cache_max) .and. &
(l <= mo_integrals_cache_max) ) then
get_mo_bielec_integral = mo_integrals_cache(i,j,k,l)
else
ii = l-mo_integrals_cache_min
ii = ior(ii, k-mo_integrals_cache_min)
ii = ior(ii, j-mo_integrals_cache_min)
ii = ior(ii, i-mo_integrals_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
get_mo_bielec_integral = dble(tmp)
else
ii = l-mo_integrals_cache_min
ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
get_mo_bielec_integral = mo_integrals_cache(ii)
endif
end
@ -390,16 +408,15 @@ double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map)
! Returns one integral <ij|kl> in the MO basis
END_DOC
integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
double precision, external :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) > mo_integrals_threshold) then
double precision, external :: get_mo_bielec_integral
!DIR$ FORCEINLINE
get_mo_bielec_integral_schwartz = get_mo_bielec_integral(i,j,k,l,map)
else
tmp = 0.d0
get_mo_bielec_integral_schwartz = 0.d0
endif
end