10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00

Introduced screening.irp.f

This commit is contained in:
Anthony Scemama 2020-05-12 18:48:51 +02:00
parent 4699ad5822
commit 08089a4dad
9 changed files with 81 additions and 50 deletions

View File

@ -0,0 +1,16 @@
logical function ao_one_e_integral_zero(i,k)
implicit none
integer, intent(in) :: i,k
ao_one_e_integral_zero = .False.
if (.not.(read_ao_one_e_integrals.or.is_periodic)) then
if (ao_overlap_abs(i,k) < ao_integrals_threshold) then
ao_one_e_integral_zero = .True.
return
endif
endif
if (ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then
ao_one_e_integral_zero = .True.
endif
end

View File

@ -85,9 +85,10 @@ double precision function get_ao_two_e_integral_erf(i,j,k,l,map) result(result)
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
integer :: ii integer :: ii
real(integral_kind) :: tmp 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 PROVIDE ao_two_e_integrals_erf_in_map ao_integrals_erf_cache ao_integrals_erf_cache_min
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then if (ao_two_e_integral_zero(i,j,k,l) < ao_integrals_threshold ) then
tmp = 0.d0 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 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 tmp = 0.d0
@ -127,10 +128,11 @@ subroutine get_ao_two_e_integrals_erf(j,k,l,sze,out_val)
integer :: i integer :: i
integer(key_kind) :: hash integer(key_kind) :: hash
double precision :: thresh double precision :: thresh
logical, external :: ao_one_e_integral_zero
PROVIDE ao_two_e_integrals_erf_in_map ao_integrals_erf_map PROVIDE ao_two_e_integrals_erf_in_map ao_integrals_erf_map
thresh = ao_integrals_threshold thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0 out_val = 0.d0
return return
endif endif
@ -156,11 +158,12 @@ subroutine get_ao_two_e_integrals_erf_non_zero(j,k,l,sze,out_val,out_val_index,n
integer :: i integer :: i
integer(key_kind) :: hash integer(key_kind) :: hash
double precision :: thresh,tmp double precision :: thresh,tmp
logical, external :: ao_one_e_integral_zero
PROVIDE ao_two_e_integrals_erf_in_map PROVIDE ao_two_e_integrals_erf_in_map
thresh = ao_integrals_threshold thresh = ao_integrals_threshold
non_zero_int = 0 non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0 out_val = 0.d0
return return
endif endif

View File

@ -291,8 +291,10 @@ subroutine compute_ao_two_e_integrals_erf(j,k,l,sze,buffer_value)
double precision :: ao_two_e_integral_erf double precision :: ao_two_e_integral_erf
integer :: i integer :: i
logical, external :: ao_one_e_integral_zero
logical, external :: ao_two_e_integral_zero
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
buffer_value = 0._integral_kind buffer_value = 0._integral_kind
return return
endif endif
@ -302,7 +304,7 @@ subroutine compute_ao_two_e_integrals_erf(j,k,l,sze,buffer_value)
endif endif
do i = 1, ao_num do i = 1, ao_num
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then if (ao_two_e_integral_zero(i,j,k,l)) then
buffer_value(i) = 0._integral_kind buffer_value(i) = 0._integral_kind
cycle cycle
endif endif
@ -618,6 +620,7 @@ subroutine compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value)
double precision :: integral, wall_0 double precision :: integral, wall_0
double precision :: thr double precision :: thr
integer :: kk, m, j1, i1 integer :: kk, m, j1, i1
logical, external :: ao_two_e_integral_zero
thr = ao_integrals_threshold thr = ao_integrals_threshold
@ -634,7 +637,7 @@ subroutine compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value)
if (i1 > j1) then if (i1 > j1) then
exit exit
endif endif
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thr) then if (ao_two_e_integral_zero(i,j,k,l)) then
cycle cycle
endif endif
if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thr ) then if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thr ) then

View File

@ -333,11 +333,10 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result)
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
integer :: ii integer :: ii
real(integral_kind) :: tmp 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_in_map ao_integrals_cache ao_integrals_cache_min
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then if (ao_two_e_integral_zero(i,j,k,l)) then
tmp = 0.d0
else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then
tmp = 0.d0 tmp = 0.d0
else else
ii = l-ao_integrals_cache_min ii = l-ao_integrals_cache_min
@ -427,9 +426,8 @@ complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result)
complex(integral_kind) :: tmp complex(integral_kind) :: tmp
PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then logical, external :: ao_two_e_integral_zero
tmp = (0.d0,0.d0) if (ao_two_e_integral_zero(i,j,k,l)) then
else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then
tmp = (0.d0,0.d0) tmp = (0.d0,0.d0)
else else
ii = l-ao_integrals_cache_min ii = l-ao_integrals_cache_min
@ -481,11 +479,10 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
integer :: i integer :: i
integer(key_kind) :: hash integer(key_kind) :: hash
double precision :: thresh logical, external :: ao_one_e_integral_zero
PROVIDE ao_two_e_integrals_in_map ao_integrals_map PROVIDE ao_two_e_integrals_in_map ao_integrals_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0 out_val = 0.d0
return return
endif endif
@ -511,11 +508,10 @@ subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val)
integer :: i integer :: i
integer(key_kind) :: hash integer(key_kind) :: hash
double precision :: thresh logical, external :: ao_one_e_integral_zero
PROVIDE ao_two_e_integrals_in_map ao_integrals_map PROVIDE ao_two_e_integrals_in_map ao_integrals_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0 out_val = 0.d0
return return
endif endif
@ -540,12 +536,13 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
integer :: i integer :: i
integer(key_kind) :: hash integer(key_kind) :: hash
double precision :: thresh,tmp 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_in_map
thresh = ao_integrals_threshold
non_zero_int = 0 non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0 out_val = 0.d0
return return
endif endif
@ -555,12 +552,12 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
integer, external :: ao_l4 integer, external :: ao_l4
double precision, external :: ao_two_e_integral double precision, external :: ao_two_e_integral
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then if (ao_two_e_integral_zero(i,j,k,l)) then
cycle cycle
endif endif
call two_e_integrals_index(i,j,k,l,hash) call two_e_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_map, hash,tmp) call map_get(ao_integrals_map, hash,tmp)
if (dabs(tmp) < thresh ) cycle if (dabs(tmp) < ao_integrals_threshold) cycle
non_zero_int = non_zero_int+1 non_zero_int = non_zero_int+1
out_val_index(non_zero_int) = i out_val_index(non_zero_int) = i
out_val(non_zero_int) = tmp out_val(non_zero_int) = tmp
@ -584,10 +581,12 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out
integer :: i,k integer :: i,k
integer(key_kind) :: hash integer(key_kind) :: hash
double precision :: tmp 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_in_map
non_zero_int = 0 non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0 out_val = 0.d0
return return
endif endif
@ -598,7 +597,7 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out
integer, external :: ao_l4 integer, external :: ao_l4
double precision, external :: ao_two_e_integral double precision, external :: ao_two_e_integral
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then if (ao_two_e_integral_zero(i,j,k,l)) then
cycle cycle
endif endif
call two_e_integrals_index(i,j,k,l,hash) call two_e_integrals_index(i,j,k,l,hash)
@ -630,10 +629,12 @@ subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,s
integer :: i,k integer :: i,k
integer(key_kind) :: hash integer(key_kind) :: hash
double precision :: tmp 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_in_map
non_zero_int = 0 non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0 out_val = 0.d0
return return
endif endif
@ -646,7 +647,7 @@ subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,s
integer, external :: ao_l4 integer, external :: ao_l4
double precision, external :: ao_two_e_integral double precision, external :: ao_two_e_integral
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then if (ao_two_e_integral_zero(i,j,k,l)) then
cycle cycle
endif endif
call two_e_integrals_index(i,j,k,l,hash) call two_e_integrals_index(i,j,k,l,hash)

View File

@ -0,0 +1,15 @@
logical function ao_two_e_integral_zero(i,j,k,l)
implicit none
integer, intent(in) :: i,j,k,l
ao_two_e_integral_zero = .False.
if (.not.(read_ao_two_e_integrals.or.is_periodic)) then
if (ao_overlap_abs(j,l)*ao_overlap_abs(i,k) < ao_integrals_threshold) then
ao_two_e_integral_zero = .True.
return
endif
endif
if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then
ao_two_e_integral_zero = .True.
endif
end

View File

@ -300,22 +300,17 @@ subroutine compute_ao_two_e_integrals(j,k,l,sze,buffer_value)
double precision :: ao_two_e_integral double precision :: ao_two_e_integral
integer :: i integer :: i
logical, external :: ao_one_e_integral_zero
logical, external :: ao_two_e_integral_zero
if (ao_overlap_abs(j,l) < thresh) then
buffer_value = 0._integral_kind if (ao_one_e_integral_zero(j,l)) then
return
endif
if (ao_two_e_integral_schwartz(j,l) < thresh ) then
buffer_value = 0._integral_kind buffer_value = 0._integral_kind
return return
endif endif
do i = 1, ao_num do i = 1, ao_num
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then if (ao_two_e_integral_zero(i,j,k,l)) then
buffer_value(i) = 0._integral_kind
cycle
endif
if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh ) then
buffer_value(i) = 0._integral_kind buffer_value(i) = 0._integral_kind
cycle cycle
endif endif
@ -1173,6 +1168,7 @@ subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
double precision :: integral, wall_0 double precision :: integral, wall_0
double precision :: thr double precision :: thr
integer :: kk, m, j1, i1 integer :: kk, m, j1, i1
logical, external :: ao_two_e_integral_zero
thr = ao_integrals_threshold thr = ao_integrals_threshold
@ -1189,10 +1185,7 @@ subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
if (i1 > j1) then if (i1 > j1) then
exit exit
endif endif
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thr) then if (ao_two_e_integral_zero(i,j,k,l)) then
cycle
endif
if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thr ) then
cycle cycle
endif endif
!DIR$ FORCEINLINE !DIR$ FORCEINLINE

View File

@ -25,7 +25,7 @@
!$OMP local_threshold)& !$OMP local_threshold)&
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
!$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, & !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, &
!$OMP ao_overlap_abs, ao_two_e_integral_alpha, ao_two_e_integral_beta) !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta)
allocate(keys(1), values(1)) allocate(keys(1), values(1))
allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), &
@ -48,8 +48,8 @@
l = ll(1) l = ll(1)
j = jj(1) j = jj(1)
if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & logical, external :: ao_two_e_integral_zero
< ao_integrals_threshold) then if (ao_two_e_integral_zero(i,k,j,l)) then
cycle cycle
endif endif
local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j)

View File

@ -28,7 +28,7 @@
!$OMP local_threshold)& !$OMP local_threshold)&
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
!$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, & !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, &
!$OMP ao_overlap_abs, ao_two_e_integral_alpha, ao_two_e_integral_beta) !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta)
allocate(keys(1), values(1)) allocate(keys(1), values(1))
allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), &
@ -51,8 +51,8 @@
l = ll(1) l = ll(1)
j = jj(1) j = jj(1)
if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & logical, external :: ao_two_e_integral_zero
< ao_integrals_threshold) then if (ao_two_e_integral_zero(i,k,j,l)) then
cycle cycle
endif endif
local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j)

View File

@ -26,7 +26,7 @@
!$OMP local_threshold)& !$OMP local_threshold)&
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
!$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, & !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, &
!$OMP ao_overlap_abs, ao_two_e_integral_alpha, ao_two_e_integral_beta) !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta)
allocate(keys(1), values(1)) allocate(keys(1), values(1))
allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), &
@ -49,8 +49,8 @@
l = ll(1) l = ll(1)
j = jj(1) j = jj(1)
if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & logical, external :: ao_two_e_integral_zero
< ao_integrals_threshold) then if (ao_two_e_integral_zero(i,k,j,l)) then
cycle cycle
endif endif
local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j)