From 08089a4dad7e83aab0f1dd6da4fce2b088751355 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 May 2020 18:48:51 +0200 Subject: [PATCH] Introduced screening.irp.f --- src/ao_one_e_ints/screening.irp.f | 16 +++++++ src/ao_two_e_erf_ints/map_integrals_erf.irp.f | 9 ++-- .../two_e_integrals_erf.irp.f | 9 ++-- src/ao_two_e_ints/map_integrals.irp.f | 43 ++++++++++--------- src/ao_two_e_ints/screening.irp.f | 15 +++++++ src/ao_two_e_ints/two_e_integrals.irp.f | 21 +++------ src/hartree_fock/fock_matrix_hf.irp.f | 6 +-- src/kohn_sham/fock_matrix_ks.irp.f | 6 +-- src/kohn_sham_rs/fock_matrix_rs_ks.irp.f | 6 +-- 9 files changed, 81 insertions(+), 50 deletions(-) create mode 100644 src/ao_one_e_ints/screening.irp.f create mode 100644 src/ao_two_e_ints/screening.irp.f diff --git a/src/ao_one_e_ints/screening.irp.f b/src/ao_one_e_ints/screening.irp.f new file mode 100644 index 00000000..aa23a9b8 --- /dev/null +++ b/src/ao_one_e_ints/screening.irp.f @@ -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 + diff --git a/src/ao_two_e_erf_ints/map_integrals_erf.irp.f b/src/ao_two_e_erf_ints/map_integrals_erf.irp.f index b3d56c41..206eea75 100644 --- a/src/ao_two_e_erf_ints/map_integrals_erf.irp.f +++ b/src/ao_two_e_erf_ints/map_integrals_erf.irp.f @@ -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 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_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 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 @@ -127,10 +128,11 @@ subroutine get_ao_two_e_integrals_erf(j,k,l,sze,out_val) 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_overlap_abs(j,l) < thresh) then + if (ao_one_e_integral_zero(j,l)) then out_val = 0.d0 return 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(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_overlap_abs(j,l) < thresh) then + if (ao_one_e_integral_zero(j,l)) then out_val = 0.d0 return endif diff --git a/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f b/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f index 97debfab..88c74ac0 100644 --- a/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f +++ b/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f @@ -291,8 +291,10 @@ subroutine compute_ao_two_e_integrals_erf(j,k,l,sze,buffer_value) double precision :: ao_two_e_integral_erf 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 return endif @@ -302,7 +304,7 @@ subroutine compute_ao_two_e_integrals_erf(j,k,l,sze,buffer_value) endif 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 @@ -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 :: thr integer :: kk, m, j1, i1 + logical, external :: ao_two_e_integral_zero 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 exit 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_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thr ) then diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 5272096d..c0ec9695 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -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 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 !DIR$ FORCEINLINE - if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then - tmp = 0.d0 - else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then + if (ao_two_e_integral_zero(i,j,k,l)) then tmp = 0.d0 else 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 PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min !DIR$ FORCEINLINE - if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then - tmp = (0.d0,0.d0) - else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then + 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 @@ -481,11 +479,10 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) integer :: i integer(key_kind) :: hash - double precision :: thresh + logical, external :: ao_one_e_integral_zero 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 return endif @@ -511,11 +508,10 @@ subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val) integer :: i integer(key_kind) :: hash - double precision :: thresh + logical, external :: ao_one_e_integral_zero 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 return 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(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 - thresh = ao_integrals_threshold 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 return 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 double precision, external :: ao_two_e_integral !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 endif call two_e_integrals_index(i,j,k,l,hash) 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 out_val_index(non_zero_int) = i 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(key_kind) :: hash double precision :: tmp + 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_overlap_abs(j,l) < thresh) then + if (ao_one_e_integral_zero(j,l)) then out_val = 0.d0 return 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 double precision, external :: ao_two_e_integral !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 endif 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(key_kind) :: hash double precision :: tmp + 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_overlap_abs(j,l) < thresh) then + if (ao_one_e_integral_zero(j,l)) then out_val = 0.d0 return 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 double precision, external :: ao_two_e_integral !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 endif call two_e_integrals_index(i,j,k,l,hash) diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f new file mode 100644 index 00000000..5095deec --- /dev/null +++ b/src/ao_two_e_ints/screening.irp.f @@ -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 diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index a2bde897..acc9ab7a 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -300,22 +300,17 @@ subroutine compute_ao_two_e_integrals(j,k,l,sze,buffer_value) double precision :: ao_two_e_integral 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 - return - endif - if (ao_two_e_integral_schwartz(j,l) < thresh ) then + + if (ao_one_e_integral_zero(j,l)) then buffer_value = 0._integral_kind return endif do i = 1, ao_num - if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) 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 + if (ao_two_e_integral_zero(i,j,k,l)) then buffer_value(i) = 0._integral_kind cycle 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 :: thr integer :: kk, m, j1, i1 + logical, external :: ao_two_e_integral_zero 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 exit endif - if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thr) then - cycle - endif - if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thr ) then + if (ao_two_e_integral_zero(i,j,k,l)) then cycle endif !DIR$ FORCEINLINE diff --git a/src/hartree_fock/fock_matrix_hf.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f index 4a62dcf3..d7d8fa7d 100644 --- a/src/hartree_fock/fock_matrix_hf.irp.f +++ b/src/hartree_fock/fock_matrix_hf.irp.f @@ -25,7 +25,7 @@ !$OMP local_threshold)& !$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_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(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & @@ -48,8 +48,8 @@ l = ll(1) j = jj(1) - if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & - < ao_integrals_threshold) then + logical, external :: ao_two_e_integral_zero + if (ao_two_e_integral_zero(i,k,j,l)) then cycle endif local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) diff --git a/src/kohn_sham/fock_matrix_ks.irp.f b/src/kohn_sham/fock_matrix_ks.irp.f index a11906db..6c3c5f9f 100644 --- a/src/kohn_sham/fock_matrix_ks.irp.f +++ b/src/kohn_sham/fock_matrix_ks.irp.f @@ -28,7 +28,7 @@ !$OMP local_threshold)& !$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_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(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & @@ -51,8 +51,8 @@ l = ll(1) j = jj(1) - if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & - < ao_integrals_threshold) then + logical, external :: ao_two_e_integral_zero + if (ao_two_e_integral_zero(i,k,j,l)) then cycle endif local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) diff --git a/src/kohn_sham_rs/fock_matrix_rs_ks.irp.f b/src/kohn_sham_rs/fock_matrix_rs_ks.irp.f index 2cfeb109..17972a79 100644 --- a/src/kohn_sham_rs/fock_matrix_rs_ks.irp.f +++ b/src/kohn_sham_rs/fock_matrix_rs_ks.irp.f @@ -26,7 +26,7 @@ !$OMP local_threshold)& !$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_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(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & @@ -49,8 +49,8 @@ l = ll(1) j = jj(1) - if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & - < ao_integrals_threshold) then + logical, external :: ao_two_e_integral_zero + if (ao_two_e_integral_zero(i,k,j,l)) then cycle endif local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j)