From 3918134a4f6151e2824f77c635933133b7d6607a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 19 Sep 2014 11:35:53 +0200 Subject: [PATCH] Added Schrwartz screening --- src/BiInts/ao_bi_integrals.irp.f | 1 + src/BiInts/map_integrals.irp.f | 17 +++++++++++++++-- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f index c4e47cf2..272f5ad8 100644 --- a/src/BiInts/ao_bi_integrals.irp.f +++ b/src/BiInts/ao_bi_integrals.irp.f @@ -215,6 +215,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax PROVIDE gauleg_t2 ao_integrals_map all_utils + PROVIDE ao_bielec_integral_schwartz integral = ao_bielec_integral(1,1,1,1) real :: map_mb diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index 04188328..a638bc39 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -62,8 +62,14 @@ double precision function get_ao_bielec_integral(i,j,k,l,map) real(integral_kind) :: tmp PROVIDE ao_bielec_integrals_in_map !DIR$ FORCEINLINE - call bielec_integrals_index(i,j,k,l,idx) - call map_get(map,idx,tmp) + 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 + call bielec_integrals_index(i,j,k,l,idx) + call map_get(map,idx,tmp) + endif get_ao_bielec_integral = tmp end @@ -92,6 +98,8 @@ subroutine get_ao_bielec_integrals(j,k,l,sze,out_val) do i=1,sze if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh ) then out_val(i) = 0.d0 + else if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then + out_val(i)=0.d0 else !DIR$ FORCEINLINE call bielec_integrals_index(i,j,k,l,hash) @@ -126,7 +134,12 @@ subroutine get_ao_bielec_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_ non_zero_int = 0 do i=1,sze + integer, external :: ao_l4 + double precision, external :: ao_bielec_integral !DIR$ FORCEINLINE + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then + cycle + endif call bielec_integrals_index(i,j,k,l,hash) call map_get(ao_integrals_map, hash,tmp) if (dabs(tmp) < thresh ) cycle